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January,  1987,  Reno,  Nevada. 
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AIAA/DGLR/JSASS  International  Electric  Propulsion  Conference,  May, 
1987,  Colorado  Spring,  Colorado. 
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4.  Shoureshi,  R.,  Pourki,  F.,  "Stability  analysis  of  electro-magneto  plasma 
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and  Control. 


Accomplishments:  Two  concurrent  studies  on  stabilizability  (Shoureshi,  Pourki)  and 

modeling  (Murthy,  Goodfellow)  of  MPD  thrusters  have  been  carried  out.  In  the  area 

of  modeling  the  main  accomplishments  have  been. 

1.  Formulation  of  a  model  for  optimized  nozzle  geometry  for  a  thermal-magnetic 
plasma  dynamic  thruster. 

2.  Formulation  of  a  model  for  energy  transfer  process  in  the  vicinity  of  electrodes  in 
a  thermal-magnetic  plasma  dynamic  thruster. 

Following  is  the  highlights  of  the  main  accomplishments  in  stabilizability  and 

observability  of  the  MPD  thrusters. 

1.  Lyapunov  stability  analysis  of  the  equilibrium  state  of  MPD  thrusters  for  the 
general  case  of  nonzero  plasma  velocity,  and  special  case  of  zero  velocity. 

2.  Stabilizability  of  the  equilibrium  state  of  MPD  thrusters. 

3.  Derivation  of  the  characteristics  of  a  general,  nonlinear  hyperbolic  system 
represented  by  partial  differential  equations  (PDE),  and  their  relations  to 
controllability  of  the  system. 

4.  Spectral  analysis  of  the  MPD  thruster  with  zero  plasma  velocity. 

5.  Study  of  observability  of  nonlinear  distributed  parameter  systems,  and  formulation 
of  observability  criteria  for  a  special  case  of  MPD  thrusters. 


MODEL  OF  PLASMA  IN  THE  VICINITY  OF  AN  ELECTRODE 


■;5 


•a 

I 

$ 

.«* 

u' 

c 

I 

I 


•» 


tjd 

iJ3 

I 

4 

•  ri 


PERSONNEL 

S.N.B.  Murthy,  Principal  Investigator 
K.  Goodfellow,  Graduate  Research  Asst. 


PUBLICATIONS  DURING  1986-1987 


1.  Pourki,  F.,  Goodfellow,  K.,  Shoureshi,  R.,  and  Murthy,  S.N.B.,  Observability  of 
Heat  Transfer  to  MPD  Thruster  Electrodes,  International  Electric  Propulsion 
Conference,  May  1987  AIAA-87-1070 

2.  Murthy,  S.N.B. ,  Shoureshi,  R.,  and  Pourki,  F.,  An  Approach  to  MPD  Engine 
Instabilities,  AIAA-87-0384. 

3.  Murthy,  S.N.B.,  Plasmadynamic  Thruster  Nozzle,  Accepted  for  1987  JANAF 
Propulsion  Meeting,  Dec.  1987. 

DISCUSSIONS  HELD 

1.  Dr.  A.  Friedman,  Dept,  of  Mathematics,  Purdue  University 

2.  Dr.  A.  Trippi,  European  Space  Agency,  Electric  Propulsion  Branch. 

3.  Dr.  R.  Hanson,  Stanford  University. 

4.  Dr.  A.  Eddy,  Georgia  Institute  of  Technology. 

5.  Dr.  T.K.  Bose,  Indian  Institute  of  Technology,  Madras. 

6.  Dr.  U.  Daybelge,  University  of  Istanbul,  Turkey. 

MAIN  ACCOMPLISHMENTS 

1.  Formulation  of  a  model  for  optimized  nozzle  geometry  for  a  thermal-magnetic 
plasma  dynamic  thruster. 

2.  Formulation  of  a  model  for  energy  transfer  processes  inthe  vicinity  of  electrodes 
in  a  thermal-magnetic  plasma  dynamic  thruster:  Details  on  this  model  are  provided  in 
the  following. 


A  MODEL  FOR  ELECTRODE  VICINITY  INTERACTIONS 


Prepared  by  K.  Goodfellow 
1.  INTRODUCTION 

In  a  number  of  problems  related  to  heat  transfer  to  electrodes,  it  has  been  usual  in 
the  past  to  neglect  the  influence  of  plasma  structure  in  the  vicinity  of  the  electrode  in 
comparison  with  the  overall  thermal  boundary  layer.  In  some  approaches  to  electrode 
surface  deterioration,  including  wear,  plasma  or  material  jets  and  diffusional  processes 
have  been  considered.  However,  most  of  those  analyses  also  neglect  the  details  of  the 
plasma  processes  in  the  electrode  region  for  which  the  characteristic  lengths  are  Debye 
length,  molecular  mean  free  path,  recombination  length  and  Larmor  radius.  S.A.  Self 
(1983,1987)  and  collaborators  have  included  consideration  of  such  near-electrode 
processes  in  the  case  of  an  isothermal  plasma  under  certain  other  approximations. 

We  have  utilized  the  S.A.  Self  approach  as  a  bench  mark  and  proceeded  to  set  up 
a  detailed  model  for  electrode  vicinity  processes  in  the  case  of  a  thermal-magnetic 
thruster  where  the  current  flow  and  the  electrode  potential  drop  and  charge  are 
obviously  of  primary  consideration. 

The  model  is  developed  in  the  context  of  energy  transfer  to  electrodes.  However, 
one  of  the  most  important  contributions  of  the  model  is  to  gain  an  understanding  of 
the  changes  in  those  processes  as  a  function  of  input  current  while  the  mass  flux  of 
heated  gas  is  held  constant.  At  some  value  of  current,  it  is  generally  accepted,  an 
instability,  called  onset  instability,  arises.  The  objective,  then,  is  to  establish  the 
nature  of  possible  measurements  in  the  electrode  vicinity  that  can  provide  an 
indication  of  the  near-onset  conditions  in  the  thruster. 

While  we  are  yet  to  quantify  the  changes,  it  is  our  belief  that  in  the  case  of  a 
cathode,  the  model  will  demonstrate  the  importance  of  (a)  the  growth  of  ion  number 
density,  (b)  the  lag  in  heating  of  ions  and  (c)  the  gradual  change  in  the  electrode 
potential. 
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2.  MODEL  DESCRIPTION 

A  number  of  zones  can  be  identified  in  the  vicinity  of  the  electrode  based  on 
characteristic  length  scales  such  as  the  recombination  length  (<r),  the  molecular  mean 
free  path  ei  n>  e  and  i  referring  to  neutral,  electron  and  ion,  respectively)  and  the 
Debye  length  (^>).  The  rate-governed  processes  in  those  zones  will  depend  further 
upon  characteristic  times.  Figure  1  presents  a  schematic  of  the  electrode  region  of  a 
MPD  engine.  A  thermal  boundary  layer  is  postulated  that  may  be  different  from  the 
momentum  boundary  layer.  The  Debye  length  and  the  mean  free  path  of  species  are 
both  small  compared  to  the  recombination  length  but  may  be  either  less  than  or  be 
comparable  to  ^  e. 

The  working  fluid  is  assumed  to  be  a  monatomic  gas  that  is  undergoing  single 
ionization  so  that  the  species  considered  are  atoms,  ions  and  electrons  obeying 
Maxwell-Boltzmann  statistics.  The  plasma  in  the  free  stream  is  assumed  to  be 
partially  ionized  with  a  possibility  of  both  thermal  and  charge  nonequilibrium.  The 
electron  and  ion  temperatures  (Te  T;)  are  therefore  assumed  to  be  different.  The  Saha 
equation  for  reaction  is  modified  in  order  to  take  into  account  non-equilibrium  and 
inelastic  collisions.  The  plasma  is  assumed  to  be  subject  to  Joule  heating  and  to 
induced  and  applied  electro-magnetic  fields. 

The  state  of  plasma  within  and  therefore  ^  may  only  be  described  in  terms  of 
collisionless  particles.  It  may  be  pointed  out  that  ^  e  and  ^  are  not  sharply  defined 
boundaries. 

Within  the  electrode,  two  regions  have  been  identified  in  Figure  1:  the  first  is  a 
region  where  melting  and  evaporation  may  be  taking  place  and  the  second,  a  region 
with  pure  conductive  heat  transfer.  The  first  is  referred  to  as  a  "mushy"  region. 


2.1.  Model  of  Plasma  in  the  Vicinity  of  an  Electrode 


A  partially  ionized  two-temperature  non-equilibrium  plasma  flow  is  considered, 
which  obeys  Maxwell-Boltzmann  statistics.  The  induced  magnetic  field  is  included  for 
consideration. 

The  plasma  is  assumed  to  be  composed  of  three  monatomic  ideal  gasses,  electron, 
ion  and  neutral.  The  pressure  of  each  gas  species  is  given  by  Ps  =  nskTs.  The  total 
gas  pressure  is  then  given  by 

P  =  £  P.  (2-1) 

9 

The  energy  of  each  species  is  the  sum  of  the  translational  and  the  internal  energies. 
The  translational  energy  of  each  species  is  given  by 

3  nskTs 


A  Cartesian  orthogonal  coordinate  system  is  used  with  coordinates  that  are  parallel  to 
the  magnetic  field,  perpendicular  to  the  magnetic  field  and  perpendicular  to  the 
magnetic  and  the  electric  fields. 
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The  electron  current  and  heat  flux  are  given  by  the  following  equations.  The 
electron  transport  properties  are  significantly  affected  by  the  magnetic  field  and 
therefore  these  effects  are  included. 


j«  =  °li  €n  +  cr±  €x  +  <%  b  X  £  +  <^iiT  V,|  Te  +  <t>x  Vj_  Te 
+  0HTtx  VTe 


(2.3) 


Z  -  -  7  —  Te  -  V„  V„  Te  -  \'±  V,  Te  -  VHt  X  VTe 


~  4>  nT  £||  —  Te  e±  —  T„  0hT  b  Xf 


(2.4) 


The  transport  coefficients  (<7,  <^T,  and  \')  are  presented,  for  example  in  Mitchner  and 
Kruger [1974],  in  the  form  of  integral  functions  of  C.  The  transport  coefficients  are  also 
available  in  other  forms  from  sources  such  as  Bose  [1987]. 

The  heavy  particle  transport  properties  are  not  affected  by  the  presence  of  a 
magnetic  field  unless  the  field  is  very  strong.  For  the  system  considered  the  magnetic 
field  is  assumed  to  be  sufficiently  weak  so  that  the  heavy  particle  transport  equations 
can  be  written  for  a  partially  ionized  plasma  without  a  magnetic  field.  Accordingly, 
the  following  relations  may  be  written  for  the  heavy  particle  transport  properties. 
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D,  =  thermal  diffusion  coefficient 


species  is  required.  Only  collisional  reactions  will  be  considered.  Specifically  the 
three-body  recombination  reaction  where  the  third  body  is  an  electron  may  be  written 
as  follows. 


Cv' 


(2.11) 


The  generation  rate  equation  is  then  given  by  the  following. 
dn. 
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=  Hinnov— Hirchberg  recombination  coefficient 


(2.12) 

(2.13) 


The  equilibrium  concentrations  are  given  by  the  Saha  equation,  namely. 
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where  g,  i?  equal  to  the  electron  energy  partition  function  and  the  ionizational 
energy.  Since  a  singly  ionized  gas  is  considered,  the  electron  and  ion  generation  rates 
are  equal  (na  =  h,). 


2.2.  Energy  Balance  in  the  Vicinity  of  a  Cathode 

As  a  particular  example  of  electrode-associated  processes,  we  examine  the  energy 
balance  in  the  vicinity  of  a  cathode.  The  equations  describing  the  state  of  the  plasma 
including  non-equilibrium  and  multi-temperature  effects  are  presented.  Referring  to 
Figure  1,  the  energy  balance  equations  applicable  to  different  regions  are  presented 
along  with  the  associated  current  flux  equation.  It  is  assumed  that  the  cathode 
material  is  catalytic  and  recombination  processes  within  the  material  give  rise  to  a 
flux  of  neutrals  into  the  Debye  region  and  also  to  heat  generation  within  the  electrode. 
A  region,  referred  to  as  the  "mushy"  region,  is  identified  in  the  electrode  material 
adjoining  the  plasma.  In  that  region  it  is  assumed  that  phase  change  processes,  solid 
to  liquid,  liquid  to  vapor  and  also  solid  to  vapor,  are  expected  to  occur. 

The  energy  equations  are  written  specifically  for  each  region  shown  in  Figure  1.  In 
all  regions  the  plasma  is  assumed  to  be  non-neutral  (na^nj)  and  two-temperature 
governed  (Ta  TJ. 

2.2.1.  Region  1  (solid) 

The  energy  equation  in  the  solid  material  is  given  by  the  well-known  Stephan  or 
heat  equation,  namely. 
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where  Eit  is  the  energy  stored  in  the  material,  (Ein  —  Eout),  net  energy  transfer 
through  the  material  by  conduction,  and  Egen,  heat  generated  within  the  material  due 
to  Joule  heating. 


The  current  flux  is  given  by 

T=  *.o.idE 


(2.16) 


The  energy  input  to  the  solid  region  is  determined  from  an  energy  balance  at  the 
surface  y  «■  0.  The  energy  balance  yields  namely, 
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where  Hm  is  the  latent  heat  of  fusion  of  the  material 


2.2.2.  Region  2  (Mushy) 

Region  2  is  the  "mushy"  region  where  solid,  liquid  and  gaseous -states  may  be 
present  simultaneously.  The  energy  for  the  "mushy"  region  is  similar  to  the  energy 
equation  for  the  solid  region  except  that  an  additional  term  is  needed  to  account  for 
the  energy  associated  with  the  material  phase  change,  Q  (T).  If  the  material  is 
assumed  to  "pure"  then  there  will  be  no  conduction  across  the  mushy  region.  The 
region  will  be  at  a  uniform  temperature  because  it  is  undergoing  a  phase  change.  If 
the  material  is  not  considered  "pure"  then  there  can  be  conduction  across  the  region 
because  there  may  be  a  small  temperature  gradient.  The  energy  transfer  in  the  mushy 
region  is  then  given  by  the  following  model  equation  [Fasano  and  Primcerio  1981], 
which  is  in  the  nature  of  a  modified  Stephan  equation. 


PmushyCPnmiAy  ^  ^  [^ushy^T  j  +  <7  ^  hy  ^  ^ 


where 

j  *“  ^mu.hy^ 

and  =  <^mushy(T) 

temperature  dependent. 


(2.18) 


(2.19) 


=  electrical  conductivity  of  the  mushy  region  which  is 


The  two  surfaces  of  the  mushy  region  at  y  =  0  and  y  =  y1(  are  not  fixed.  The 
surface  at  y  =  0  is  allowed  to  recess  into  the  solid  material  as  melting  occurs.  Melting 
will  occur  in  the  region  as  long  as  the  energy  input  to  the  region  (at  the  surface  y  —  yt 
and  from  Joule  heating)  exceeds  the  energy  used  in  the  phase  transition  plus  the 
energy  removed  to  the  solid  region  by  conduction.  If  the  energy  inputs  and  outputs 
are  equal  the  boundary  will  not  recede.  The  surface  at  y  =  yl  will  recede  because  of 
the  evaporation  of  material  at  the  surface. 

The  energy  input  to  the  mushy  region  from  the  plasma  is  determined  from  an 
energy  balance  evaluated  at  the  surface  y  =  yj,  as  shown  in  Figure  1.  The  cathode  is 
assumed  to  have  a  catalytic  surface  where  incident  ions  and  electrons  recombine  and 
are  re-emitted  as  neutrals.  Thermionic  emission  may  also  occur  if  the  cathode 


temperature  is  sufficiently  high.  It  should  be  noted  that  the  thermionically  emitted 
electrons  provide  an  additional  localized  current  which  in  turn  produces  an  additional 
localized  Joule  heating.  The  localized  heating  may  lead  to  localized  evaporation  or  to 

the  eruption  of  material.  The  overall  energy  balance  equation  is  then  as  follows. 

* 

Q  “b  Qcond  Qp  d“  q^d  d~  Qgur  (2.20) 

where  qp  is  the  net  energy  transfer  to  the  surface  from  the  plasma  due  to  the  kinetic 
energy  and  potential  energy  of  the  particles,  qj.ad,  the  net  radiation  to  the  surface  from 
the  plasma,  Qsut,  the  energy  generation  at  the  surface  due  to  recombination  of  ions 
and  electrons  Q  ,  the  energy  associated  with  phase  change  and,  q^^,  energy 
conducted  into  the  mushy  region. 

2.2.3.  Region  3  (Sheath) 

The  region  immediately  adjacent  to  the  electrode  and  within  a  distance  of  the 
order  of  Debye  length  (^>)  of  the  surface  is  the  sheath  region.  In  this  region,  charge 
separation  occurs  and  a  net  negative  charge  exists  because  of  an  excess  number  of 
electrons.  The  region  is  considered  collisionless  in  the  sense  that  only  electron-neutral 
and  ion-neutral  collisions  are  present.  These  collisions  are  included  because  of  the 
large  number  of  neutrals  in  the  region.  Since  neutrals  being  emitted  from  the  surface 
do  not  experience  a  force  from  the  fields,  they  tend  to  remain  near  the  surface.  They 
are  removed  through  diffusion  driven  by  the  concentration  gradient.  All  other 
collisions  are  assumed  to  be  negligible. 

The  energy  transfer  to  the  surface  from  the  plasma  is  from  the  collision  of  particles 
on  the  surface.  Since  a  cathode  is  being  considered  the  particles  of  interest  are  the 
ions.  The  energy  of  the  particles  is  in  two  forms,  the  kinetic  energy  due  to  their 
motion  and  the  potential  energy  associated  with  moving  charged  particles  through  an 
electric  potential.  The  electric  potential  will  tend  to  move  the  electrons  away  from  the 
electrode  while  accelerating  the  ions  towards  the  electrode. 

Since  the  interest  is  in  particles  that  strike  the  surface,  the  velocity  component 
normal  to  the  surface  is  the  one  of  interest.  The  describing  energy  equations  are  given 
by  the  following.  (Mitchner  and  Kruger  1974] 


The  species  heat-flux  vector,  "qs,  is 

q,  -  JX  7  n,  m,  C2  C±  f,  d3c  (2.21) 

00  4 

and  the  particle  potential  energy,  Ws,  is 

W*  =  lyl y,  (2-22) 

where  the  electric  potential  Os  is  given  by  the  following  Poisson's  equation. 

=  /_*  f,f,d3c  (2.23) 

The  species  current-flux  vector,  j3,  given  by 

T,  =  J_l  n.  C1  C1  f3  d3c  (2.24) 

The  distribution  fs  is  determined  from  the  Boltzmann  equation,  namely 
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dfs  <9u~ 

~  C/jDs  ac7  ^7  = 


where  is  the  rate  of  increase  of  the  property  of  interest  (mass,  momentum,  energy 
or  charge)  due  to  collisions  between  particles  of  type  s  with  particles  of  type  k. 


2.2.4.  Region  4 


Region  4  is  similar  to  region  3,  as  they  are  both  considered  collisionless.  The 
describing  equations  are  therefore  the  same  for  the  two  regions.  Since  this  region  is 
outside  of  the  sheath,  it  is  expected  to  contain  a  greater  number  of  ions  than  region  3. 
The  electron-neutral  and  ion-neutral  collisions  also  become  negligible  because  of  a 
large  decrease  in  the  number  of  neutrals  in  the  region.  The  collision  parameter  in 
equation  2.25  is  therefore  equal  to  zero. 


2.2.5.  Regions  5  and  6 


Regions  5  and  6  are  assumed  to  be  collision-dominated  and  can  therefore  be 
described  using  the  hydrodynamic  approximation.  Particle  recombination  is  assumed 
to  absent  in  Region  5  because  it  is  within  a  distance  of  the  order  of  the  recombination 
length  of  the  electrode  surface.  The  recombination  is  included  in  Region  6.  The 
resulting  energy  balance  equations  for  Region  5  are  as  follows. 
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The  energy  equations  for  Region  6  are  given  by  the  following.  The  recombination 
energy  is  included  in  the  electron  energy  equation. 
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2.2.7.  Region  7  (Free  Stream) 


In  the  free  stream  region,  the  various  gradients  are  assumed  to  very  small 
compared  to  the  other  terms  and  are  neglected.  The  rate  of  change  of  energy  within 
the  system  is  equal  to  the  generated  energy  from  Joule  heating  minus  the  energy  lost 
from  radiation,  recombination  and  collisions.  The  energy  equations  can  therefore  be 
simplified  as  follows. 


d  3  ,  _  2me  _  3  ( 

dt  ne  I  2  kTe  +  £i  I]  ~  Je‘^e  ~  mh  '/ehIle  I  k  lT 


(Te  -  Th)  -  R 


a  f  nhkT"  ^  f  k(Tt-T.)-R, 


Je  =  en  +  ei  +  °h  b  X  e 


(2.30) 


(2.31) 


(2.32) 


3.  PROPOSED  SCHEME  FOR  UTILIZATION  OF  THE  MODEL 

The  vicinity  of  an  electrode,  in  particular  the  cathode,  in  an  MPD  thruster  is 
considered  according  to  the  simplified  model  presented  in  Figure  2.  The  geometry  and 
the  flow  rate  of  propellant  are  considered  fixed.  The  thruster  is  then  expected  to  be 
operated  at  two  current  levels,  namely  a  low  current  level  in  which  the  plasma  is 
partially  ionized  and  a  high  current  level  at  which  the  plasma  is  nearly  fully  ionized. 

The  objective  is  to  establish  a  methodology  for  determining  the  plasma  properties 
in  the  cathode  vicinity  and  for  relating  changes  in  such  properties  with  respect  to 
current  to  the  occurrence  of  onset  conditions. 
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The  describing  equations  for  the  model  are  considered  in  two  forms,  a  collision- 
dominated  plasma  and  a  collisionless  plasma.  The  sheath  and  free-fall  zones  are 
considered  as  collisionless,  since  they  fall  within  a  distance  less  than  one  mean  free 
path  (4  e)  from  the  electrode  surface.  Plasma  that  is  located  at  y  >  fie  is  considered 
as  collision  dominated  plasma.  For  the  purposes  here  Regions  3,  4  and  5  from  Figure  1 
are  considered,  that  is  the  plasma  within  the  recombination  length  of  the  cathode 
surface. 


Assumptions: 

1.  One  dimensional  model  (normal  to  the  cathode  surface). 

2.  Maxwell-Boltzmann  statistics  apply  to  the  particles. 

3.  Isotropic  properties. 

4.  Negligible  convection. 

5.  Negligible  magnetic  field  (induced  and/or  applied). 

6.  Cathode  absorbs  all  incident  species  (ions  are  reemitted  as  neutrals). 

7.  Monoenergtic  particles  (ions  and  thermionic  electrons). 

8.  Thermionic  electrons  are  emitted  at  the  cathode  temperature. 

9.  No  recombination  or  ionization  occurs  in  y  <  ^. 

10.  The  energy  of  the  thermionic  (beam)  electrons  is  negligibly  small  in  the  collision 
dominated  region. 

11.  Collisions  are  included  for  the  region  £i  e  <  y  <  €e -r 

12.  Diffusion  neglected  in  the  y  <  ^  e. 

3.1.  Govering  Equations 


3.1.1.  Collisionless  Region 

The  sheath  and  free-fall  regions  are  considered  to  contain  a  collisionless  plasma, 
that  is,  no  interparticle  collisions  occur  within  these  regions.  The  total  energy  of  each 
particle  is  therefore  constant  since  energy  cannot  be  gained  or  lost  by  a  particle 
without  a  collision.  Three  species  of  particles  are  considered,  positive  ions,  plasma 
electrons  and  thermionic  (or  beam)  electrons.  Neutral  particles  will  also  be  present. 

The  solution  procedure  is  to  solve  for  the  number  density  of  each  species  using  the 
energy  and  continuity  equations  and  then  use  Poissons  equation  to  solve  for  the 
electric  field  (E)  and  the  electric  potential  ((f)).  Once  the  electric  potential  is 
determined,  the  number  densities  and  and  velocities  of  each  species  can  be  determined. 

Electrons  produced  by  thermionic  emission  at  the  cathode  surface  are  accelerated 
through  the  region  by  the  electric  field.  Their  total  energy  is  composed  of  two  parts, 
the  particle  kinetic  energy  and  the  potential  energy  in  the  form: 


Eb  =  —  mbvb2  —  e<f>  —  constant  (3.1) 

2 

The  beam  electrons  are  assumed  to  be  emitted  with  temperatures  equal  to  the  cathode 
temperature  (Tc).  The  initial  beam  electron  energy  is  therefore: 

Ec  =  J-W  =  |kT»  (3.2) 

Using  equations  3.1  and  3.2,  the  beam  electron  energy  can  be  written  in  the  form: 

ymbvb2  =  e(0c  -  <t>)  +  Ec  ,  (3.3) 

which  shows  that  the  kinetic  energy  is  equal  to  the  change  in  potential  energy  plus  the 
initial  kinetic  energy.  The  initial  kinetic  energy  is  expected  to  be  very  small  compared 
to  the  other  terms  but  is  retained  to  prevent  a  possible  singularity  from  occuring  at 
the  electrode  surface  ( 4>  =  d>c).  The  potentials  are  used  as  absolute  values 
(<t>  =  \—<j\  ,  d>c  ==  I  — 0CI ).  The  flux  of  beam  electrons  within  the  region  is  continuous, 
since  no  ionization  or  recombination  occurs.  The  continuity  equation  yields: 

jb  =  enbvb  (3.4) 

Solving  equations  3.3  and  3.4  for  the  number  density,  nb  yields: 


Jb  mb  (.  A  „ 
nb  -  — '  —  ’  e  4>c-  4>  +  Ec 

C  At  > 


The  thermionic  emission  current  is  assumed  to  be  dependent  only  on  the  cathode 
temperature  and  material  and  is  given  by  the  empirical  relationship: 

_  5040  0 

Jb  -  (PA)  T2  10  T  (3.6) 

where  pA  and  0  are  material  properties  and  T  is  the  electrode  temperature  in  degrees 
Kelvin  [Smithells  1952].  For  the  purposes  here,  the  effect  of  space  charge  on 
thermionic  emission  has  been  neglected,  that  is,  a  space  charge  in  the  vicinity  of  the 
cathode  acting  either  to  enhance  or  to  inhibit  the  electron  emission  is  neglected. 


The  ions  are  assumed  to  be  monoenergetic  with  an  initial  kinetic  energy  of  — mjVjo2 

2 


at  the  plane  where  they  enter  the  collisionless  region  (y  =  4e)-  The  ion  energy 
equation  can  be  written  as: 

ymjVi2  -  e<t>  =  Eio  (3.7) 

where  Eio  is  the  initial  ion  energy  at  y  =  fj  e 

Eio  =  ~"mivio2  “  e<^2  (3.8) 


Ei0  =  e^o  (3-9) 

where  <t>0  is  the  equivalent  potential  drop  for  the  ion  to  obtain  the  energy  Eio.  From 
equations  3.7  and  3.9  the  ion  velocity  can  be  determined  as: 
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The  ion  distribution  is  continuous  due  to  the  absence  of  collisions  and  continuity 
yields: 


J«  =  en^i  =  eniovic 


(3.11) 


where  nio  and  vio  represent  the  ion  number  density  and  velocity  at  the  plane  where  the 
ions  enter  the  collisionless  region  at  y  =  ^  e.  Solving  for  n;  yields: 
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The  plasma  electrons  are  assumed  to  be  Maxwellian  particles.  The  number  density 
is  given  by: 

—  ed>  +  E, 

(3.13) 


ne  =  ne2exp 


kTe 


where  E,*,  is  the  electron  kinetic  energy  at  the  free-fall  edge  y  =  Cx  e,  the  potential 
energy  is  ed>2,  and  the  number  density  is  ne2.  The  energies  are  related  such  that: 


Eeo  -  e02  =  0 

Equation  3.13  can  be  rewritten  as 
-  e  -  <£2  j 


(3.14) 


ne  =  ne2  exp 


kT. 


(3.15) 


The  electron  current  in  the  collisionless  region  is  considered  in  two  parts.  The  electrons 
that  move  towards  the  cathode  with  initial  energies,  Eeo  <  e(^>c  —  (f>2)  are  repelled  and 
return  to  the  collision  dominated  region.  The  electrons  with  Eeo  >  e(<j!>c  —  02)  cross 
the  collisionless  region  and  strike  the  cathode.  These  high  energy  electrons  make  up 
the  electron  current  je. 
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(3.16) 


The  solution  of  Poisson’s  equation  in  the  sheath  region  will  be  different  from  the 
solution  in  the  free-fall  region  because  of  charge  separation  within  the  sheath. 
Poisson’s  equation  is  used  to  solve  for  the  electric  potential  distribution  (0). 
Substituting  the  particle  number  densities  into  Poisson’s  equation  yields: 
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.  .  .  the  guessed  value  of  <*<  are  known  constants,  equation  3.18  can  be 

Since  j„  Jb  and  the  g  d<i>  ,  th  lven  boundary  conditions 

integrated  once  using  the  integrating  facto,  —  and  the  gtven 

0  «  <?3  and  E  =  Ej  at  y  -  <b  “  y'*M: 
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to  write  equation  3.19  in  the  form: 


Poisson’s  equation  can  be  used  again  to  wnte  equation  — 
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Equation  3.20  can  then  be  numerically  integrated  to  determine  the  electric  potential 
distribution  within  the  sheath.  Once  the  potential  is  determined,  the  number  densities 
of  each  species,  the  electric  field  and  the  space  charge  can  be  determined  using 
equations  3.5,  3.12,  3.13,  and  3.17  respectfully. 


A  separate  equation  is  needed  to  determine  if  the  guessed  value  of  <f>c  is  correct. 
This  equation  is  obtained  by  performing  a  force  balance  on  the  ions.  This  yields 
equation  3.24  which  has  a  resemblance  to  Ohms  law  for  a  plasma  with  collisions. 
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which  can  be  simplified  for  the  case  at  hand  where  j*  —  constant,  to  yield 
mj  jj2  dn, 
e3  Qi3  dy 


(3.25) 


The  derivation  of  equation  is  given  in  appendix  I.  A  new  value  for  </>c  is  obtained  by 
appling  equations  3.20  and  3.25  at  the  cathode  surface.  This  procedure  is  repeated, 
starting  again  from  y  ==  ^e,  until  the  guessed  and  calculated  values  of  4>c  are  within  a 
chosen  tolerance. 


The  solution  of  Poisson’s  equation  is  simplified  within  the  free-fall  region  due  to  the 
absence  of  charge  separation,  that  is  pe  being  constant.  The  value  of  pc  is  determined 
by  the  charge  at  the  edge  boundary  (y  =  ^  e).  Poisson’s  equation  can  be  integrated 
once  to  yield: 


_iiL_E_E  _£L[(  _v| 

dy  E-Es  e„  T'  yl 


dy 


and  integrated  again  to  yield: 
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3.1.2.  Collision  Dominated  Region 


The  governing  equations  for  the  collision  dominated  region  are  given  by  the 
hydrodynamic  plasma  equations.  These  equations  are  applicable  when  the  number  of 
collisions  within  the  plasma  are  sufficiently  large  so  that  average  properties  can  be 
assigned  to  the  plasma  properties.  The  equations  are  further  simplified  here  for  the 
one  dimensional  case.  The  derivation  of  these  equations  can  be  found  in  a  number  of 
references  such  as  Mitchner  and  Kruger[l974]. 
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(3.29) 


j.  =  enjUj 

(3.30) 

j«  -  eneUe 

(3.31) 

Jb  "  enbub 

(3.32) 
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Poisson’s  equation 

(3.33) 
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These  equations  are  used  to  solve  for  the  variables  ne,  ni(  nn,  nb,  Te,  Th,  j„  j;,  p  ,  Ui( 
U„  Un,  Ub,  E  and  <t>. 


The  characteristic  length  scales  are  given  by  the  equations: 


n«Q.e 


(3.42) 


where  Qie  is  the  collisional  cross-section.  The  sheath  thickness  is  determined  using  the 
Child-Langmuir  law  [Chen  1985). 
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Alternatively,  or  the  following  equation  can  be  used. 
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The  sheath  thickness  is  equal  to  the  larger  of  the  two  values  for  Cq  from  equations  3.43 
and  3.45. 
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4.  PLANNED  SOLUTION  PROCEDURE 

The  overall  solution  procedure  is  outlined  in  the  flow  chart  in  Figure  3.  The 
general  solution  is  found  by  using  the  solutions  for  each  of  the  regions.  The  governing 
equations  are  solved  for  each  region  starting  with  the  collision-dominated  region  and 
proceeding  towards  the  cathode  surface.  The  parameter  values  at  the  boundaries  are 
determined  by  the  preceding  region  and  are  used  as  the  beginning  boundary  conditions 
for  the  next  region.  The  process  is  repeated  until  the  estimated  values  for  collisionless 
region  potential  drop  and  beam  electron  number  density  agree  with  the  calculated 
values.  Once  convergence  is  achieved,  the  distributions  of  the  various  parameters  can 
be  studied. 

Four  test  cases  are  expected  to  be  examined.  They  are  for  combinations  of  high 
and  low  current  density  levels  and  high  and  low  cathode  temperatures.  It  is  expected 
that  estimation  of  distributions  of  the  critical  parameters  can  be  determined  by 
examining  the  results  of  these  four  cases. 
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B 

c* 

0 

t 

r 

T 

q 

if 

0 

cp 

e 

E 

E 

h 

j 

k 

m 

M 

n 

P 

* 

Q 

R 
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Description 
magnetic  field 

particle  velocity  in  laboratory  reference  frame 

peculiar  particle  velocity 

applied  electric  field 

body  force  on  particle 

current  flux 

heat  flux 

mass  velocity 

diffusion  velocity  in  fluid  reference  frame 
specific  heat 
electron  charge 
energy 

electric  field  (l-D) 

Planck  constant 
current  density  (l-D) 

Boltzmann  constant 
mass 

gas  molecular  weight 
number  density 
pressure 

phase  change  energy 
radiation  energy 
temperature 


E 


V 

velocity  (1-D) 

W 

potential  energy 

X 

position 

r 

general  electric  field 

«i 

ionization  potential 

permittivity  constant 

n i 

viscosity 

X,  V 

thermal  conductivity 

0 

electric  potential 

0T 

thermal  diffusion 

p 

density 

<7 

electrical  conductivity 

P 

mobility 

mean  collision  frequency  of  type  s  particles  with  type  k  particles 

T 

shear  stress 

& 

particle  charge 

4 

Debye  length 

e 

mean  free  path 

Subscripts 

b 

thermionic  (beam)  electrons 

c 

cathode 

cond 

conduction 

e 

electron 

:«ii 
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k 

mushy 


JL 

H 
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Definitions 


b  -2. 
B 


species  of  type  k 
mushy  region 
neutral 

reference  state 
species  of  type  s 
radiation 

solid  material  region 
surface 

parallel  to  magnetic  field 
perpendicular  to  magnetic  field 

mutually  perpendicular  to  magnetic  and  electric  fields 

directions 

free  stream  value 

boundary  y  =  ^ 

boundary  y  =  €x  e 

boundary  y  = 


C,=  mean  particle  speed  =  j - 

[  ^ 

?  =  ^  +  uXB 
-  AP 

r  =  E*  + — - 

nee 

f,  =  velocity  distribution  function 
ps  =  n3m,  =  species  density 


mmmmmrmsm 

-  23  - 


APPENDIX  I 


A  relationship  between  current  and  electric  field  for  a  collisionless  region  can  be 
developed  as  follows.  The  total  current  is  given  by 


j  =  Ji  +  Jb  —  Je 


(1.1) 


where  the  electron  and  beam  currents  are  constant  for  a  given  cathode  temperature 
and  collisionless  region  electric  potential  drop. 


jb  =  (PA)  T2  10 


S040 

T 


Je  = 


Aeo 


'  8kTe  ' 

'k 

■  exp 

^  J 

~  e(&  ~  ^2) 
kT„ 


(1.2) 


(1.3) 


ji  =  enjVj 

A  force  balance  on  the  ion  yields 


(1.4) 


dVj  dv, 

F  =  —  eE  =  m— —  =  mv:  — — 
dt  1  dy 


Substituting  for  v;  yields 


(1.5) 


(1.6) 


Figure  1. 


Planned  Solution  Procedure 


Solution  procedure. 


Figure  3.  Continued 


Calculate  Change  in  Collision- Dominated 
Region  From  Beam  Electrons 


Modify 

Collision- Dominated 
Region  for  Beam 
Electron  Influence 


STABILIZ  A BILITY  ANALYSIS  OF  MFD  THRUSTERS 


1.  INTRODUCTION 

One  of  the  main  concerns  of  the  control  theory  for  systems  governed  by  partial 
differential  equations,  is  whether  or  not  the  system  is  stable,  if  unstable,  how  the 
system  could  be  stabilized  by  applying  a  proper  set  of  control  inputs  to  the  system. 
These  systems  which  are  often  called  distributed  parameter  systems,  are  represented 
by  a  set  of  states  which  are  functions  of  time  and  spatial  coordinate(s).  At  each 
instance  of  time,  every  state  of  the  system  has  a  distribution  on  the  spatial  plane  (i.e. 
the  state  belongs  to  an  infinite-dimension  functional  space).  However,  for  finite 
dimensional  systems  governed  by  a  set  of  ordinary  differential  equations,  a  state  is 
defined  by  a  scaler  at  every  instance  of  time.  For  finite  dimensional  systems  Lyapunov 
stability  theorem  lias  become  an  important  vehicle  in  derivation  of  stability  analysis  of 
dynamics  of  the  system.  This  approach  attempts  to  make  statements  about  stability 
of  motion  of  a  dynamic  system  without  any  knowledge  of  the  solutions  to  its  governing 
equations. 

Although  the  development  of  Lyapunov’s  stability  theorem  for  ordinary  differential 
equations  has  been  widely  investigated,  its  application  to  solutions  of  partial 
differential  equations  (distributed  parameter  systems)  has  been  limited.  Many  stability 
results  for  distributed  parameter  systems  have  been  derived  by  use  of  approximation 
methods.  These  methods,  in  general,  use  reduction  of  the  partial  differential  equations 
by  discretization  |1|*  or  by  a  truncation  of  the  modal  expansion  ;2,3j  to  a  finite  large 
order  set  of  ordinary  differential  equations.  Such  approaches  may  not  give  sufficient 


conditions  for  system  stability.  Another  advantage  of  Lyapunov’s  method  over 
approximate  methods  is  that  Lyapunov's  stability  theorem  can  be  applied  to  both 
linear  and  nonlinear  systems. 

The  attempt  to  apply  Lyapunov's  method  to  partial  differential  equation  (PDK) 
was  made  by  Massera  and  Zubov  .4,5..  Massera  derived  sufficient  conditions  for 
stability  of  equilibrium  solutions  (steady  states)  of  system  of  PDE.  The  generalization 
of  Lyapunov's  stability  theorem  based  on  the  existence  of  a  Lyapunov  functional  was 
established  by  Zubov.  He  derived  necessary  and  sufficient  condition  for  the  stability  of 
invariant  set  of  dynamical  systems  in  general  metric  spaces.  The  application  of 
Lyapunov  stability  theorem  based  on  the  work  of  Zubov  has  been  investigated  by 
several  authors.  Hsu  6  applied  this  theory  to  a  nuclear  reactor  system.  Wang  [ 7 j 
considered  the  stability  of  those  evolution  equations  whose  solutions  involve  a 
semigroup  property.  There  are  also  many  other  applications  which  utilize  Lyapunov 
functions  directly  to  study  special  problems  (8,91-  A  completely  rigorous  and  abstract 
approach  to  the  theory  of  Lyapunov  stability  for  infinite  dimensional  system  was 
studied  by  10  and  11  .  Stability  study  of  an  equilibrium  solution  of  a  magneto¬ 
plasma-dynamic  (MPD)  system  for  the  special  case,  where  the  plasma  equilibrium 
velocity  is  zero  was  addressed  by  the  authors  12.  During  the  past  budget  year,  the 
study  has  been  focused  on  the  problem  of  controlling  the  general  equilibrium  state  of 
MPD  thrusters.  Moreover,  the  requirements  for  stability  of  a  general  equilibrium  state 
of  the  system  with  an  initial  perturbation  is  derived.  This  analysis  provides  the 
foundation  for  development  of  required  control  inputs  which  guarantee  stability  of  the 


-Argument  of  denotes  the  reference  number. 


-  32  - 


MPD  system  perturbed  from  a  general  equilibrium  state.  The  problem  of 
stabilizability  which  is  tied  to  the  controllability  property  of  a  distributed  parameter 
system  has  been  studied  in  general  by  [  1 3|  and  jl4j  in  Hilbert  and  Banach  spaces, 
respectively.  Applications  of  these  studies  on  specific  systems  are  investigated  by  (loj 
and  16  for  wave  equations  with  boundary  and  internal  control  problems.  The 

generalized  approach  to  the  problem  is  investigated  and  the  results  are  presented  in 
the  following  sections. 

Section  2  presents  some  minimum  mathematical  preliminary  background  related  to 
the  treatment  of  the  problem.  The  focal  point  in  this  section  is  to  give  the  notion  of 
equivalent  norm  and  a  theoretical  basis  of  why  the  bilinear  functional  (i.e.  Lyapunov 
functional)  can  be  regarded  as  an  equivalent  norm.  In  section  3  the  semigroup 
properties  of  solutions  of  systems  of  PDE's  is  studied,  while  in  section  4  the 

significance  and  the  relation  of  these  properties  to  stability  conditions  is  presented.  In 
section  5  the  simple  models  of  MPD  thrusters,  as  it  is  studied  in  [ 1 7 j  and  {18},  is 

derived.  The  control  oT  equilibrium  states  of  the  MPD  system,  represented  by  the 

models  to  avoid  singularity  at  choking  and  stability  conditions  for  the  convergence  of 
perturbation  to  the  equilibrium  state  have  been  studied  and  presented  in  sections  6 
and  7,  respectively. 

The  report  is  concluded  by  presenting  two  possible  approaches,  based  on  the  theory 
of  characteristics  and  ^election  of  some  energy  function,  which  could  lead  to  the 
determination  of  general  system  stabilizability  requirements.  These  approaches  can  be 
regarded  as  basis  for  future  research  in  applying  the  stabilizability  and  controllability 
f  o  t  he  MPD  system. 
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2.  PRELIMINARIES 


This  section  provides  the  mathematical  basis  and  preliminaries  for  the  following 
sections.  The  objective  is  to  provide  the  readers  with  a  quick  review  of  the  required 
definitions  and  theorems. 

Definition:  Let  X  be  a  vector  space  over  real  or  complex  field  F.  A  norm  on  X, 

denoted  by  II  •  II  is  a  real-valued  function  on  X  with  the  following  properties: 

a)  II  x  II  >  0  if  x  =*=  0  and  II  x  II  =  0  for  x  =  0  for  all  x  e  X 

b)  II  or  x  II  —  |  or  |  II  x  II  .  for  or  c  F  (F  =  '3i  or  '€  ) 

c)  II  x  4-  y  II  <  II  x  II  -F  II  y  II  .  for  x  and  y  e  X 

Norms  can  be  constructed  in  different  ways. 

If 

x  =  Jxj  j  i=l,2...n. | 

Then  for  infinite  dimensional  space  X,  n  — ►  x  the  following  is  defined. 


-  :u  - 


Based  on  the  definition  used  for  the  norm,  different  types  of  normed  linear  spares 
ran  be  specified.  The  two  most  commonly  used  normed  spaces  are  Hilbert  and  Banach 
spaces.  Hilbert  space  is  a  normed  linear  space  with  the  norm  being  defined  as  II- 1 lLi  or 
INI,-.  Respectively,  Banach  spaces  are  defined  with  II- 11^  or  INI,  norms. 

Lyapunov  stability  theory  evaluates  the  stability  with  respect  to  a  norm.  This 
theory  is  concerned  primarily  with  dissipative  property  of  operators,  which  can  be 
related  to  inner  product  property  of  the  space  of  state  variables.  Therefore  the 
natural  space  for  stability  analysis  would  be  space  with  inner  product  norm  (i.e. 
Hilbert  space).  However,  for  the  study  of  stability  in  general  Banach  space,  norm  can 
be  defined  in  term  of  semi-inner  product.  Since  the  stability  properties  are  invariant 
under  equivalent  norms,  the  related  definition  and  theorems  are  given  in  the  following. 

Definition:  Let  H,  with  norm  INI,2  =  <  .,.  >,  and  H,  with  norm  IHI22  =  <  ...  >2  be 
Hilbert  spaces  consisting  of  the  elements  of  a  linear  vector  space  H  and  the  given 
norms.  The  inner  products  are  called  equivalent  if  and  only  if  their  induced  norms  are 
equal.  Hence  by  this  notion  of  equivalent  norm,  as  will  be  seen  later,  one  can  preserve 
the  stability  properties  from  one  Hilbert  space  into  another  Hilbert  space  with 
equivalent  norms. 

Linear  F unctional:  Every  functional  T  :  x  — is  called  a  linear  functional  provided: 

Xj,  x.j  (  X  C  If  (Banach  Space) 

T(x2)  +  T(x2)  =  T(x,  +  x,) 

T(orx)  =  <*T(x) 

T(0)  =  0 

A  linear  functional  T  is  bounded  if 


Riesz  Representation  Theorem: 


Every  bounded  linear  functional  in  a  Hilbert  space  H  can  be  written  in  the  form 
<x,y0>  where  y0eH  is  uniquely  determined  from  T(x). 

T(x)  =  <  x,y0  >  V/xeH 

Lax-Milgram  Theorem: 

Let  H  be  a  Hilbert  space  and  let  B(x,y)  be  a  complex  valued  functional  defined  on 
the  product  Hilbert  space  H  >  H  which  satisfies  the  following  conditions: 

i.  Sesqui-linearity ,  i.e., 

B  |  (oqx,  +  e*2x2),y  |  =  a,B(x,,y)  +  a2B(x2,y) 

and 

>3  (x.  ((llyl  +  P2y2)  |  =  j3{  B(x,  y ,)  +  JJ2  0(xt,y2) 

where  /?.,  are  complex  conjugates  of  respectively, 

ii.  Boundedness,  i.e.,  there  exists  a  positive  constant  7  such  that 

|B(x,y)|  <  7llxlldlyll 

iii.  Positivity,  i.e.,  there  exists  a  positive  constant  6  such  that 

B(x,x)  >  <51  lx II2 

Then  there  exists  a  uniquely  determined  bounded  linear  operator  SeL(H,H)  with  a 
bounded  linear  inverse  S  '(LOLII)  such  that 

<x.y  >  —  B(x,Sy);  1ISII  <  1  / 

and 


,y);  US  ’ll  <  7 

This  theorem  can  be  used  in  relation  with  equivalent  inner  products  and  norms. 
Theorem:  Two  inner  products  defined  on  a  real  linear  vector  space  H  are  equivalent  if 
and  only  if  there  exists  a  symmetric  bounded  positive  definite  linear  operator  SeL(H.H) 
such  that 


<x.y>2  =  <x,sy>j 

where  the  indices  indentify  the  inner  product  in  Hilbert  spaces  1  and  2.  From  the  fact 
that  B(x,y)  =  <x,  Sy>(  =  <x,v>2,  is  sesqui-linear,  bounded  and  positive  due  to  the 
properties  of  S,  it  follows  that  <x.y>2  satisfies  all  the  properties  of  an  inner  product. 
Hence  the  norms  defined  by  these  inner  products  are  equivalent.  Application  of  the 
equivalent  norms,  yields  derivation  of  the  stability  properties  of  the  general  operators 
(A)  in  differential  equations  of  the  form 

x  =  Ax , 

from  the  knowledge  of  the  properties  of  the  operator  A.  The  concept  of  equivalent 
inner  products  can  be  extended  to  complex  Hilbert  spaces. 
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3.  SEMI-GROUPS  AND  ABSTRACT  EVALUATION  EQUATIONS 


Finite  order  systems  can  be  described  by  a  set  of  ordinary  differential  equations. 
The  simplest  systems  are  linear  autonornus  systems  which  can  be  formulated  by 

v  =  Av  v(0)  =  v0  (3-1) 

where  v:[0,Tj  — ►  is  the  solution  of  the  system  which  is  a  vector  of  continuously 

differentiable  functions.  The  matrix  A  is  a  linear  operator  which  acts  on  vct/t”  (n- 

dimensional  vector  space),  or  Ac  f  (.'/?“).  For  v  =  eAlv0  is  the  homogeneous 

solution  which  satisfies 

77(eAtvo)  =  AeAtv0  =  A(v) 

(It 

where 

a  2  4  n 

e  A  =  I  +  A  -f  T  _ T  4"  •  •  • 

2!  A! 

Usually  external  effects  such  as  body  forces  can  be  modeled  into  the  above  mentioned 
system  by  a  vector  valued  function  g  :  tf  j(),T)  — * 

v  =  Av  +  g(t)  v(0)  =  v0  (3-2) 

Then  the  solution  is 

t 

v(t)  =  eAtv0  +  /  eA(t"s>  g(s)ds  (3-3) 

O 

Usually  the  excitation  term  g(t)  can  be  represented  in  terms  of  a  vector  valued  (input) 
function  u(t)  :  jO.Tj  — ►  Jim  with  g(t)  =  Bu(t)  where  B  :  ■/fn  — *  .y?n 

Abstract  Equations: 

In  the  case  of  distributed  parameter  systems,  the  mathematical  description  is 
usually  given  by  a  set  of  partial  differential  equations.  In  order  to  generalize  the 
results  from  the  finite-dimensional  systems,  the  set,  of  partial  differential  equations  can 


be  transformed  into  abstract  form  of  equations  (3-1)  and  (3-2).  The  state  of  system  v 
belongs  to  some  Banach  space  B  or  more  commonly,  to  a  Hilbert  space  H.  In  this  case 
the  abstract  equations  represent  a  set  of  infinite-dimensional  states  and  the  operator  A 
can  no  longer  be  treated  as  a  bounded  operator.  Because  it  is  a  differential  operator. 
For  finite-dimensional  systems  the  following  properties  are  considered. 

(a)  v(0;vo)  =  v0  (i.e.  state  at  time  0,  starting  with  initial  data  v0) 

(b)  v(tt  +  t2;  v0)  =  v[tp  v  ( 1 2 ;  v„) )  =  v  1 12;  v(t(;  vQ)  ] 

(c)  v(t,v0)  is  continuous  in  t  and  v0  (i.e.  if  one  of  them  changes  slightly,  then  the 
solution  should  not  be  changed  drastically). 

Based  on  abstraction  of  the  above  conditions,  the  following  definition  can  be  given. 

Definition: 

A  (strongly  continuous)  semigroup  is  an  operator  T(t):  R+  — *  (JE  (V),  where 
J'  (V)  is  the  space  of  all  bounded  operators  on  V  into  V.  The  following  properties 
characterize  a  semigroup 

(a)  T(o)  =  1 

(b)  T(t,  +  t2)  =  T(t,)T(t2) 

(c)  iim  T(t)v  =  v,  for  all  '■<  V  (i.e.  'f  is  stronglv  continuous  at  t,  =  0). 
t— o 

Theorem: 


Let  T ( t )  be  a  semigroup  on  H*  to  7  ( \ ' ) ,  where  V  is  a  Banach  Space,  then 


-  :w  - 


(a)  ilT(t) II  is  bounded  on  every  compact  interval  of  0,  cxj!  such  that  IIT(t)ll  <  Me"'  for 
some  M,  c 0€'Jt  . 

(b)  T(t)  is  strongly  continuous,  on  0,  ooj. 

Definition:  The  infinitesimal  generator  of  a  strongly  continuous  semigroup  T(t)  is 

defined  by  A.  when 

Av  =  lirn  I  — •  (T(t)v— v) 
t— o  (  t 

v  e  D(A)  C  V 

It  should  be  noted  that  in  general,  A.  will  be  an  unbounded  operator.  The  operator  A 
is  closed  (i.e.  its  range  and  domain  coverage  to  some  element  of  their  respective 
spaces),  and  its  domain  covers  the  closure  of  space  V  (i.e.  it  is  dense  in  V). 

Theorem:  Let  T(t)  be  a  strongly  continuous  semigroup  on  a  Banach  space  V  with 
infinitesimal  generator  A.  If  v„  <  D(A),  then 

(a)  T(t )  v„  c  D( A)  for  t  >0 

(b)  ~(T(t)vn)  =  A(T(t)  vG)  =  T(t)  Av0,  t>() 
dt 

(V)  ~(T(t)  vj  -  Ac(T(t)  vj  -  T(t)An  v„  for  v()  r  l>(An);  t  >0 


(d)  T(t)v0  -  v0  =-  /  T(s)Av„  ds;  t>() 

O 

(e)  D(An)  is  dense  in  V  for  n  1.  2,  ...  and  A  is  closed. 


B 


For  finite  dimensional  system,  Laplace  transform  of  eAI  is  L(cAI)  —  (SI— A)  '.  This  can 
be  generalized  to  semigroups  by  the  following  proposition 
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Proposition:  IF  T(t)  is  a  strongly  continuous  semigroup  with  infinitesimal  generator  A, 
then  Re(s)  >  w  for  se  p( A) 

where,  p(A)  =  resolvent  set  of  A  is  defined  by 


p[\)  =  js;(SI— A)  1  fJt'  (V);  bounded  linear  operator  onV  f. 


Note  that  ui  is  given  by 


w  =  Inf  WIIT(t)ll  <  Me-'1;  M,  f. 


Hille- Yoshida  Theorem: 


(a)  A  is  a  closed  linear  operator  on  V  such  that  D(A)  is  dense  in  V. 


(b)  (XI— A)  1  exists  for  \  >uj  where  X,t cf.Ji  . 


(r)  ll(Xl-A)-,nll  <  —  -  ;  X>w  :  m  -  1,  2 . 

X — u/tl 

then  A  generates  a  strongly  continuous  semigroup  T(t)  with  the  norm  I  IT  ( t )  1 1  <  Me"'1. 


Definition:  Consider  strongly  continuous  semigroup  T(t)  with 

IIT(t)ll  <  Me"'1  M  >0,  u)  <  oc 

if  c^=0  then  I  IT  ( t )  1 1  <  M  and  T(t)  is  called  an  equi-bounded  semi-group  with  t  >  0. 
Definition:  For  equi-bounded  semigroup  when  M  =  1, 

IIT(t)ll  <  1,  t  >  0  , 

then  T(t)  is  called  contraction  semigroup. 

In  the  stability  theory  of  semigroups  the  contraction  semigroups  are  very 
important.  These  contraction  semigroups  are  closely  related  to  the  dissipative 
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property  of  infinitesimal  generator  A  of  the  semigroup  T.  To  study  dissipative 


property  of  operator  A  one  needs  to  apply  it  to  inner  products,  hence  Hilbert  space 


would  be  a  natural  space  for  the  study  of  dissipative  property  of  operators. 


Definition:  An  operator  A  defined  on  D(A)  C  H  (Hilbert  Space)  is  dissipative  if 


Re<Av,  v)  <  0  for  every  vd)(A] 


Phillips  and  Lumer  Theorem 


Let  A  be  a  linear  operator  with  domain  D(A)  and  range  of  A  [Ran  (A)]  both  in  the 


Hilbert  Space  H  where  D(A)  =  H  (i.e.  domain  of  A  denses  in  H),  then  A  generates  a 


contraction  semigroup  on  H  if  and  only  if  A  is  dissipative  with  respect  to  the  inner 


product  defined  on  H,  also  Ran(l-A)  -  H. 


The  proof  is  avoided  as  in  earlier  theorems,  however,  it  should  be  mentioned  that  if 


the  hypothesis  of  the  above  theorem  holds,  then  one  can  derive 


II - - - II  <  X" 1 

(XI— A)  - 


for  all  X  >  0,  hence  the  Hille- Yoshida  theorem  could  be  used.  Namely,  by  setting  M 


—  1.  ur=0  one  can  find 


1  IT  ( t )  I !  <  Me1*'1  =  1 


Which  is  the  fact  that  the  operator  A  generates  a  contraction  semigroup. 


This  result  can  be  generalized  to  Banach  Spaces  if  a  Semi-inner  product  is  used  for 


inner  product. 


( 'orollary:  If  A  is  a  closed  linear  operator  with  dense  domain  in  H,  then  A  generates  a 


•ontraction  semigroup  if  and  only  if  A  and  A*  are  dissipative.  Note:  A*  is  the  adjoint 


operator  of  A  given  by: 


<x.Ay>  =  <A*x,y> 

So  far  the  theory  has  covered  the  properties  of  an  operator  A  with  regard  to 
semigroups.  In  general  evolution  equation,  operator  A  can  represent  a  hyperolic  or 
parabolic  problem.  However  the  specific  problem  which  will  be  considered  here  is  a 
hyperbolic  system  of  the  class 

v  =  Lv  ve(L2(n),En) 

The  states  form  an  Euchidean  vector  space  (i.e.  v  is  a  vector  valued  function) 
where  each  element  of  the  vector  forms  a  function  which  belongs  to  infinite 
dimensional  Hilbert  space  defined  on  one  dimensional  spatial  region  (domain) 
fl  =  xe[0/  j.  The  operator  L  is  similar  to  the  operator  A  discussed  before,  and  is 
unbounded  with  specific  form: 

Dv 

Lv  =  A  ~  +  Bv 
<vx 

It  is  possible  to  establish  the  solution  of  this  system  of  evolution  equations  in  terms 
of  semigroups  using  a  theorem  to  be  seen  later.  The  solutions  of  this  general  system, 
depending  on  the  conditions  specified  by  the  theorem,  can  result  in  either  groups  or 


semigroups. 


4.  STABILITY  ANALYSIS  FOR  DISTRIBUTED  PARAMETER  SYSTEMS 


In  general,  stability  theory  is  an  attempt  to  make  some  assessments  about  the 
motion  of  a  system  without  apriori  knowledge  of  the  solutions  to  its  governing  dynamic 
equation.  This  phenomena  which  is  an  intrinsic  property  of  a  system  in  a  heuristic 
sense  is  related  to  boundedness  of  the  motion  of  the  system  about  some  prescribed 
regions.  In  order  to  deline  stability  property  for  a  generalized  dynamical  system,  it  is 
appropriate  to  give  definitions  of  such  systems  and  their  invariant  sets. 

Definition:  Let  C  be  a  closed  subset  of  a  complete  metric  space.  A  dynamical  system 
on  C  is  a  family  of  maps  S  such  that  the  following  conditions  hold: 

(a)  For  any  finite  value  of  vjO),  the  vector  function  v(t,v(0))  is  defined  for  all 
t  t  (—00,  Too)  and  v(0,  v( 0 ) )  =  %(()). 

(b)  The  vector  f unction  v(t,  v(0))  is  continuous  in  the  aggregate  of  its  arguments. 

(c)  For  all  values  of  t  and  r; 

v ( t  +  r.  v(0))  =  v(t.  v(r,  v(0))) 

Hence  it  can  be  seen  that  in  general  a  dynamical  system  forms  a  non'inear  group  S 
on  C.  For  a  dynamical  system  on  C,  the  positive  semi-orbit  through  v  is  the  set 

7(v)  =  | S(t)v:t  >  0  . 

This  set  defines  the  trajectory  of  solution  with  respect  to  initial  condition  v.  The 
metric  on  the  space  U  can  lie  assigned  as  the  norm  between  two  states  v|l;  and  v|t  as: 

'*( v|ta.  v|,t)  =  [|  v|,2  -  v|Jj 

where  v|t  for  a  distributed  parameter  system  is  a  set  of  real  valued  functions  at  time  t. 


V 


-  It  - 


or 


when  d  v|tj,v|ti 


v,(t,  x)  [  i = 1 . n  and  x  =  spatial  coordinate 


vlt  =  S(0  vo 

=  0  the  state  functions  are  identical 


Definition:  An  equilibrium  state  v  of  a  dynamical  system  is  an  element  of  the  state 
space  C  such  that  d(S(t)v  ,  veq)  =  0  for  all  t  >  0.  The  set  of  all  equilibrium  states 
will  be  called  the  equilibrium  set  f  . 


Definition:  An  invariant  set  of  dynamical  system  is  a  subset  K  C  C  such  that  for  any 
initial  state  v0  e  K.  the  trajectory  v  =  S(t)v0  for  all  t  >  0  remains  in  the  set  K,  i.e.,  if 


v0  e  K,  (v  =  S(t)  v0)  c  Iv.  For  example 


'eq 


and  6  are  invariant  sets  of  the  system. 


The  distance  (metric)  of  a  state  v  from  an  invariant  set  K  is  defined  by 


d(v.K)  =  inf  d(v.v') 

v'.  K 

Definition:  An  invariant  set  K  of  a  dynamical  system  is  stable  with  metrics  d,  if  for 
every  e  >  0  there  exists  a  6(t,  t0)  such  that  if  d(vD  K)  <  <)(c,  t0)  then  d(S(t)  vQ,  K)  <  e 
for  all  t  >  tQ. 


The  invariant  set  is  asymptotically  stable  if  d(S(t)vn,  K)  — ►  0  as  t  — *  -Tex) 


Stability  Theorem  for  Distributed  Parameter  System 

Zubov  Theorem  I:  A  necessary  and  sufficient  condition  for  an  invariant  set  K  of  a 
distributed  parameter  dynamical  system  v  defined  on  space  C  to  be  stable  is  that 
there  exists  a  real  functional  L(v)  having  the  following  properties: 


p  Jf* 
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(a)  L(v)  is  defined  for  all  t  j;  0  and  all  in  some  region  with  d(v,K)  <  r. 

(b)  For  each  sufficiently  small  real  r/,  >  0  there  exists  a  real  r/2  >  0  such  that 
L(v)  >  r/2  for  d(v,K)  >  r/lt  i.e.  L(v)  is  positive  definite  with  respect  to  the  norm 


(c)  —  L(v)  <  0  for  t  >  0,  i.e.  L(v(t))  is  non-increasing  with  respect  to  all  t  >  tH,  for 
dt 

d(v0.  K)  <  (V0,  where  <*)0  is  a  sufficiently  small  positive  number. 


Zubov's  Theorem  (11):  An  invariant  set  of  K  of  a  dynamic  system  on  space  C  is 
asymptotically  stable  if  and  only  if  there  exists  a  real  functional  L(v)  having  the 
properties  (a),  (b).  (c)  of  theorem  I  and, 


(d)  —  L(v)  -—*0  as  !  — ►  +;x>.  provided  that  d(vot  K)  <  b0.  where  b0  is  a  sufficiently 
small  real  positive  number. 


Remark:  In  the  Zubov  Theorem,  the  assumption  that  the  solutions  v  satisfy  the 
properties  of  a  dynamical  system  is  needed  to  prove  the  necessary  part  of  the  theorem. 
For  the  general  system,  the  existence  of  Lyapunov  functional  provides  the  sufficiency 
condit  ion. 


Stability  of  Homogeneous  Linear  Systems 

Consider  the  abstract  linear  evolution  equation  of  the  previous  section; 

v  =  Av  v(0)  —  v„  (4-1) 

where  the  operation  A  generates  a  semigroup  T(t)  on  Banach  Space  V.  For  this 

system,  one  can  define  a  stronger  criterion  for  the  stability  by  defining  the  notion  of 

exponential  stability. 


'iv'C  '.’C-.S''. -LI*' 


Definition:  The  system  with  the  above  evolution  equation  (4-1)  is  exponentially  stable 
if  there  exists  real  positive  M  and  u,' such  that 

||T(t)||  <  Me  vt,  t  ^  0. 

Note  that  exponential  stability  implies  asymptotic  stability  but  the  converse  is  not 
always  true. 

Clearly  a  linear  system  is  exponentially  stable  if  and  only  if  its  spectrum  cr(A)  is 
contained  in  the  left  half  of  complex  plane  such  that 

Sup  Re  of  A)  <  0 

where  spectrum  is  defined. 

rr(A)  =  X:  (XI— A)v  =  0.  which  implies  (Xl— A)  is  not  one  one 

Hence  —  ui  <  0  can  be  chosen  as  sup  Re  cr(.\)  and  ||T(t)||  <  Me  ‘x't. 

Because  some  operators  do  not  have  eigenvalues  then  spectral  analysis  of  operators 
in  determination  of  their  spectrum  is  not  always  achievable.  In  these  cases  the 
Lyapunov  method  can  be  adopted  to  analyse  the  stability  of  the  solutions. 

Lyapunov  Method  for  Stability  of  Linear  Systems 

Let  operator  A  in  the  evolution  equation  (4-1)  be  generator  of  a  semigroup  T(t), 
then  the  null  solution  of  (1-1)  is  asymptotically  stable  if  there  exists  a  Lyapunov 
functional  L(v)  such  that  L(v)  ■  0  and  L(v)  <  —  7  ||v||2  for  vtD(A). 

If  the  ||‘||  is  a  Hilbert  norm  with  D(A)  =  11.  then  any  Lyapunov  functional  with 
bilinear  form  L:  (II  ■  II)  — *  ft.  by  Lax-milgram  theorem  can  be  made  into  an  inner 
product  of  t lie  form:  L  :  (v  •  w)  — *  -  v,  Sw>,  where  S  is  symmetric  real  bounded 


operator. 


’•'I 


Hence  an  equivalent  norm  can  he  defined  | |v  1 12  =  <v,  v>2  =  <v,  Sv  >  and 


tt<v,  v> 


<  I lv  I |o  <  /V<v,  v>, 


By  this  equivalence  relation,  l. lie  existence  of  Lyapunov  functional  would  lead  into 
existence  of  operator  S  with  aforementioned  properties.  Then  the  resulting  Lyapunov 
functional  by  hypothesis  satisfies; 


L(v)  =  L(v,v)  =  <v,  Sv> 


L(v)  =  lim  L(T(t)v,  T(t)v)  -  L(v.v)[-  l/t 


=  lim  |  L((T(t )  +  l)v,  (T(t )  -  I)v  |  •  l/t  =  2L(v,Av)  =  2<v,SAv>  =  2<v,Av>2 


L  II..II2 


L(v)  =  2<v.Av>,  v.  — 7||v||"  V  ~  — 


\<v,v>2  —  <  Av.v>2  =  •  (XI  -  A)v.v  >2  <  1 1(  XI  —  A  )v  2  |  v  1 12 


<V.V>.  X  -  ||(XI-A)||,  ||v||| 


where  from  the  hypothesis 


IIXI  -X 


Av.v 


X  +  -2- 


|(  XI  —  A ) 


I  rorn  this  relation  if  the  use  of  llille- Yoshida  theorem  is  made,  it  can  be  shown  that 
operator  A  is  the  generator  of  a  semigroup  T'( t )  such  that 
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which  is  a  condition  for  exponential  svmptotic  stability. 
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Recent  increase  in  space  missions  and  construction  of  the  space  station  has 
attracted  attention  to  new  alternatives  to  chemical  propulsion  systems.  One  such 
system  is  categorized  and  known  as  electrical  propulsion  thruster.  In  general, 
electrical  rockets  should  be  able  to  develop  considerably  higher  specific  impulses  than 
chemical  or  nuclear  ones.  However,  this  gain  in  specific  impulse  requires  massive 
energy  conversion  mechanism.  To  avoid  this,  the  electrical  rockets  generally  have  low 
thrust  for  navigation  in  low  gravitational  holds. 

The  propellant  of  an  electrical  rocket  consists  of  either  charged  particles  which  are 
accelerated  by  electrostatic  forces,  or  a  stream  of  electrically  conducting  fluid  (plasma) 
which  is  accelerated  by  electromagnetic  and/or  pressure  forces. 
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The  electrical  accelerators  are  divided  into  three  main  categories:  electrostatic, 
electrothermal  and  electromagnetic  accelerators. 

Electromagnetic  accelerators  use  the  conductivity  of  the  plasma  propellant  to 
create  an  electromagnetic  accelerating  force  as  a  "body  force"  within  the  plasma.  The 
main  focus  of  this  work  is  aimed  at  modeling  and  analysis  of  steady  magnetogas- 
dynamic  flow  accelerators,  among  other  categories  of  electromagnetic  accelerators. 
Figure  (5-1)  shows  a  schematic  of  such  a  system.  The  flow  of  ionized  gas  enters  the 
thruster  and  is  subjected  to  an  electric  field  E  and  a  magnetic  Held  11,  which  are 
perpendicular  to  each  other  and  to  the  gas  velocity.  Electromagnetic  acceleration 
process  is  an  aggregate  of  effects  from  compressible  gas  dynamics,  ionizer)  gas  physics, 
electromagnetic  Held  theory  and  particle  electrodynamics.  The  individual  analytic 
complexity  of  each  of  these  phenomenon  adds  to  the  level  of  difficulty  in  adequate 
theoretical  model  for  this  composite  system.  Analytical  progress  normally  stems  from 
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simplified  models  which  preserve  the  essential  physical  aspects  of  the  specific  situation. 


The  description  of  the  motion  of  plasma  in  terms  of  Maxwell-Boltzmann 
distribution  function  is  too  detailed  to  be  useful  for  many  practical  problems  in  the 
electromagnetic  acceleration  process.  In  these  cases  the  ionized  gas  medium  can  be 
considered  as  a  continuum  fluid  whose  macroscopic  physical  properties  may  be 
described  by  conservation  laws  and  Maxwell  equations.  These  governing  equations  for 
the  motion  of  plasma  will  be  gas  dynamic  equations  with  the  interaction  terms  due  to 
electromagnetic  forces.  With  this  approach  a  simplified  model  for  the  problem, 
depicted  schematically  in  figure  (5-1),  can  be  derived.  As  shown  in  this  figure  the 
plasma  is  flowing  through  the  constant  area  channel  along  the  x-coordinate.  The 
channel  is  formed  by  two  opposite  conducting  walls  connected  to  cathode  and  anode 
poles  respectively.  Between  these  walls  an  electric  field  E(x)  is  maintained  in  a  y- 
direction.  Normal  to  this  electric  field  is  a  magnetic  field  JB(x),  applied  in  the  z- 
direction.  The  model  is  considered  to  be  one  dimensional  i.e.  only  variations  in  x- 
direction  are  considered.  The  bulk  properties  of  the  gas  (plasma)  is  shared  by  all 
species  contained  in  the  plasma.  Application  of  the  conservation  laws  to  the  flow 
through  an  element  of  volume  results  in  tfie  following  equations.  From  the  continuity 
equation,  one  can  find 


<>P  +  <Hpn 

<h  cTx 

where  p  and  u  ar(.'  the  density  and  the  velocity  of  gas  respectively.  The  momentum 
equation  results  in: 


L  =  0 


(5-1) 


<Hp")  ,  <>(P" ') 

i><  +  '  iht 


ihi 

ih 
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where  p  is  the  sum  of  thermodyi.  mic  pressure  and  radiation  pressure  1’^,  however,  l’R 
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is  assumed  negligible.  Fe  is  electromagnetic  force  (Lorentz  force)  per  unit  volume,  i.e. 
Fe  =  j*B  and  Fv  represents  combination  of  viscous  and  collisional  force  and  assumed 
to  be  negligible  compared  to  other  terms.  After  expansion  of  the  momentum  equation 
and  substitution  from  continuity  equation,  the  following  is  resulted. 

p—  +  p u—  =  -  —  +  jB  .  (5-2) 

r/t  r/x  (hx. 

Due  to  the  one-dimensional  assumptions,  the  variations  in  the  y  and  z  directions  are 
negligible,  and  the  velocity  components  in  y  and  z  directions  would  be  deleted.  From 
the  energy  equation  for  the  element  of  volume  one  can  derive  the  following: 

+  itoti  _  jE  +  21hI  +  ifi. 

<7  t  (7X  <7X  (JX 


where 


specific  total  energy  (internal  energy  +  potential  energy 
-  kinetic  energy) 

specific  total  enthalpy  (h  =  e  +  P/p ) 


— -  shear  stress  tensor 


O  —  heat  llux  due  to  convection  and  radiation. 


On  the  right,  hand  side  of  the  above  energy  equation  the  term  jE  represents  the 
Joule  heating  due  to  the  application  of  electric  field.  This  term  can  be  considered  as 
the  dominant  form  of  dissipation  of  energy.  The  plasma  can  be  assumed  as  a  perfect 


gas,  with  the  state  equation 


P  -  pin’,  R  =  RA/m 


where  temperature  is  denoted  by  T,  and  RA  and  m  are  the  gas  constant  of  plasma  and 
molecular  weight  of  the  plasma,  respectively.  As  a  result,  the  specific  energy,  e,  and 
specific  enthalpy,  h,  can  be  represented  as 


wmmm 


when  the  potential  energy  terms  are  negligible.  The  specihc  heat  coefficients  cv  and  cf 
can  be  considered 

<\-  =  3/2  R 

cp  =  5/2 R 

Hence  the  energy  equation  can  be  approximated  and  re  written  as  follows: 


<)  ir  ,  <)  „•» 

-  /’C,r  +  ,.T  +-^  +,T  =  jE 

The  Maxwell  equations  can  be  written  as 


V  x  Bvec  =  Mej  +  /ie  — 


V  <E  =  -  ~ 

()t 

V  •  B  =  0 


v  ■  E  =  -P0 


For  many  practical  problems  one  can  show  that  the  displacement  current  and 

n  t 

excess  electric  charge  pp  are  negligible  and  (hat  the  energy  in  the  electric  field  is  much 
smaller  than  that  in  the  magnetic  field.  As  a  result  the  set  of  Maxwell  equations  will 
become 


mm 
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V  •  B  =  j 


v  x  e  -= 


<h 


(S-i) 


V  •  B  —  0 


V  •  E  -  0 


If  one  specifies  that  the  magnetic  field  generated  by  the  current  Rowing  in  the  gas 
is  negligible  with  respect  to  the  applied  field  B(x),  then  the  electromagnetic  field 
relations  (set  of  equations  5-4)  can  be  decoupled  from  the  dynamics  of  the  plasma. 

By  Ohm’s  law,  one  can  relate  current  density  j  with  the  applied  fields; 

j  =  o(E  -  uB),  o  =  o{p, T) 

where  rr  is  a  transport  coefficient  and  it  is  called  the  electrical  conductivity. 

After  substitution  for  in  the  energy  equation  (5-3),  one  can  arrive  at  the 

<7t  at 

following  set  of  equations  (5-5  to  5-7)  as  the  set  of  dynamical  governing  equations  for 
the  MPD  thruster 


<)P  ,  V(pu) 
()t  <h 


i)u  (hi 

~rr  +  u 
<dt  ()x 


=  0 


!i  <HpT)  +  jB 
n  i  )x  (> 


in  irr 

<)\  <• 


■  ,, /rr  i  J(E  ~  llB) 
dx  <rx  f>  cv 


(5-5) 

(5-6) 

(5-7) 


where  j  is  given  by  Ohm’s  law. 


i 


6.  CONTROL  OF  MI’D  THRUSTER  AT  EQUILIBRIUM  STATE 

In  the  set  of  non-steadv  equations  derived  in  the  previous  section,  i.e.  equations  (5- 
5  to  5-7),  it  can  be  seen  that  the  parameters  B  and  E  can  be  regarded  as  input 
functions  to  the  system.  One  of  the  crucial  questions  about  the  behavior  of  this 
system  is  how  to  characterize  the  relation  between  the  system  response  and  those 
inputs.  The  general  response  characterization  of  these  systems  which  change  with 
both  time  t  and  spatial  coordinate  x,  is  a  complex  problem.  However,  one  can  break 
the  problem  into  steps,  by  first  trying  to  make  assessments  about  equilibrium  states  of 
the  system  of  the  partial  differential  equations.  Hence,  one  can  pose  a  question  of 
what  choice  of  inputs  would  steer  the  system  to  an  equilibrium  set  of  states,  and  under 
what  conditions  the  steering  would  not  be  plausible.  In  order  to  answer  these 
arguments  about  the  equilibrium  states  of  the  system  one  has  to  derive  the  equilibrium 
set  of  equations  from  the  original  partial  dilferential  equations.  If  states  of  the  system 
are  represented  by  a  vector  valued  function  v,  as  follows 

P(t.x) 

V  =  ll(t.x) 

T(t.x). 

then  the  equilibrium  vector  would  be  denoted  by  v0 

/MM 

up(x) 

T„(x) 

Similarly  the  input s  (control)  functions  can  be  defined  by  a  vector  valued  func  tion 

U, 

l 
1 


U  - 


JB 

JK 


At  equilibrium,  system  states  reach  their  steady  state  values  and  their  time 

(  )y 

derivatives  would  be  zero,  i.e.  — 1  =0-  ;“>d  equations  (5-5  to  5-7)  will  turn  to  the 

dt 


following  equations. 


,|  dpe  pP  due 

dx  °  °'  dx  ufi  dx 


d  Up  |t  d(peTe)  lT! 

uP  - ;  +  .  — 

dx  p ,.  dx  pe 


RI’  du„ 


dTe  U2  -  upU, 


Substituting  equation  (6-1)  into  (6-2)  would  result  in  the  following  equations, 


RTe  du0  dTe  U, 

-  —  +  R  —  — 

u.  dx  dx  pe 


RTp  ditp  dTe  U2  -  ueU, 

—  +  ’Jo-]  — 

cv  dx  dx  Pe  cv 


d"p  ^  / ueU i  -  (t-1)U2  ^6_4j 

,lx  Pel^e  -  7RTe) 

d  l'e  — Ugl't  +  (iip2  —  Rljl'z  ^  -  ■ 

d\  Peuecv(Up  -  7Rd'p) 

where  7  is  the  specific  heat  ratio  cp/cv.  The  speed  of  sound  ap  is  defined  bj 
a2  _  q.RTp.  (6-4)  and  (6-5)  can  be  rewritten  into  a  vector  differential  equation  form. 


7up  7-1 


n;  -  RTP 

P0CV  fW'v 


Therefore,  for  any  equilibrium  state  configuration,  it  is  possible  to  arrive  at  a 
control  vector  in  the  local  sense  with  respect  to  the  coordinate  x.  However,  in 
transition  from  subsonic  How  to  supersonic  How,  where  ue  =  ae  (i.e.  Mach  number  is 
unity),  the  control  inputs  U,  and  1'2  become  dependent  on  each  other,  or  one  should 
choose  Ut  and  U2  according  to  a  certain  relation  such  that  the  passage  from  subsonic 
to  supersonic  velocities  and  vise  versa  would  be  plausible. 


7ueUi 

u2  =  7 - 77-  at  ue  =  7RTe 

(7-1) 


E  =  nrue 

B  “  (7-1) 

This  effect  is  the  same  as  choking  condition  in  gas  dynamics  which  is  extended  to 
magneto  gas  dynamics  ,1ft’. 


Based  on  the  definition  of  Mach  number  M  as  M  =  - - ,  it  can  be  shown  that 


VI,  «.■  ,  dM2  duo 

M*  = -  and  - —  =  — r 

7BTp  M2  u; 


dM-  _  , 

M2  ”  2  iiP 


1  clTc 


M2  fix  u„  fix  T„  fix 


if  the  substitution  for  - -  and  - -  from  equation  (6-1)  and  (6-5)  into  the  above 

dx  dx 

equation  is  made,  one  could  derive  the  following  results 


d\r  i\\-  +  1  , . 

-  ' Cl 


2  L  — L\f/2  +  1 

O  ' 


M2  dx  peuecpTe  "  pe 

For  subsonic  flow  M  <  1,  increase  in  M  is  possible  provided  the  following  inequality 


is  satisfied. 


/VipcpTe 


h\1Z  +  lj 


[^M2  +  lj  (7-0 


For  supersonic  flow  M  >  1.  increase  in  M  along  the  thruster  occurs  if, 


[7M2  +  lj 


(6-10) 


At  sonic  condition  the  ratio  of  control  inputs  -  should  comply  with  the 


aforementioned  quantity. 


Clearly,  for  the  decelerating  flow  (when  M 


decreases  along  the  flow)  the  direction  of  inequalities  (6-9)  and  (6-10)  would  be 
reversed. 

One  approach  to  construct  acceptable  control  inputs  U,  and  U2,  which  is  based  on 
exclusion  of  singularity  at  choking,  is  to  consider 


(7-1  )M2  +2j-  7»o 

^2  —  ^  \  '  •>  /  ,  v  '•■It 

17M“  +  U  (7-0 

where  is  required  to  have  the  following  properties  for  accelerating  flow; 


(6-11) 


I'  ,  >  0  lor  M  <  1 

L*3(M)  =  i:,  =  0  for  M  =  1 

l1.,  <  0  fo<-  M  >  1 


(6-12) 


Substitution  for  l7.?  from  equation  (6-11)  and  (6-12)  in  equations  (6-4)  and  (6-5)  would 


result  in  the  following 


du„  _  'Jeui  (4~1)  hi3 

dx  pe(7M2-|  l)  q,pe(M2-l) 


dTe  _  -‘2Ut  (M2— 1  /qf)  U3 

dx  />PR(7iM24 1)  pecvue(M2  — 1) 


In  order  to  avoid  singular  (unacceptable)  control  distributions,  once  should  consider  U3 
a  function  with  the  following  form; 

L'3  =  (I-M2)U4  (6-13) 

with  U,  >  0  for  accelerating  flow,  the  requirements  of  (6-12)  will  be  satisfied. 

Substitution  of  (6-13)  into  (6-1)  and  (6-5)  will  result  in  the  following  equations. 


dx  pe(7M2+l) 


dx  prR(7Mz  +  l) 


(4-1)  U4 


(M2— 1/4) 


PPev  u„ 


(6-14) 


(6-15) 


Since  the  Joule  heating  (U2)  is  a  positive  function,  for  supersonic  flow,  U3  becomes 
negative  and  from  condition  (6-1)  one  would  have; 


(M-— 1)(4'M2+1) 


for  M2  >  1 


(6-16) 


<  1,  U2  in  equation  6-1  1  is  positive  for  any  U,  >  0.  Also,  the  current  density  j 


should  be  positive  (in  the  chosen  direction  to  accelerate  the  flow),  hence  from  Ohm’s 


law  one  can  find  the  following  conditions. 


.1  -  a{E  -  uB)  >  0 


E  >  uB 


(6-17) 


l  ':>  7'  nel.'j 

From  equations  (6-11)  and  (6-13),  can  he  substituted  into  inequality  (6-17),  hence; 


for  M  <  1 


T>p^i 


(7iM2-l) 


(6-18) 


for  M  >  1 


l;i  -  lh-l)M2 


- 2 -  <  _  -Li-/  _  - 1-1  -  (6-19) 

7'iel-h  (7M2  +  1)(M2  -  l)  M2  -  1) 


It  is  clear  that  in  subsonic  flow,  inequality  (6-18)  is  always  satisfied,  hence  there  is  no 
constraint  on  U4  in  this  regime.  However,  for  supersonic  flow,  inequality  (6-19) 
represents  a  more  restrictive  constraint  on  U4  than  inequality  (6-16).  Hence  a  proper 
choice  of  U4  can  be  selected  to  satisfy  inequality  (6-19).  The  steady  state  response  can 
be  found  from  equations  (6-11)  and  (6-15)  for  some  arbitrary  choice  of  U,  >0  and  by 
selection  of  U,  according  to  the  aforementioned  process.  It  would  be  an  interesting 
proposition  to  apply  the  theory  of  optimal  control  to  liiul  the  optimal  control  inputs  U, 
and  U4  among  an  arbitrary  class  of  functions.  Those  optimal  Ult  IJ4  result  in  some 
minimized  performance  index  or  cost  function.  As  an  example,  the  cost  function  can 
be  selected  from  the  group  of  "fuel  optimal''  problems,  i.e.. 


X  < 

II'  =  /  [l'f  F  Ko,,^,)2  ^ 


Having  in  mind  that  the  optimal  control  problem  is  involved  with  control  constraint  of 


the  form  of  inequality  (6-19)  while  f  ,.  I",  an  positive  quantities. 


Perturbation  of  nonlinear  nmi-st  rn  dy  ('qua  I  ions 


Based  on  the  equilibrium  state  (represented  by  ve(x)),  it  would  be  interesting  to 
evaluate  the  system  behavior  in  the  neighborhood  of  its  equilibrium  states.  If  the 
perturbation  of  states  with  respect  to  equilibrium  is  denoted  byv,  then, 

v  —  v  +  vp  = 

P  «  Pe 
u  «  Up 
f  <K  Tp 

By  substitution  of  v  into  the  dynamic  equations  (5-5)  to  (5-7),  one  would  find; 


+  (PU  e  +  “Pe  +  fV‘e  +  P»)  =  0 

fit 

riti  ,  ,,  fi(Up+u)  R  fi(peTe  + />TC  +  pj  +  pT)  U, 

— —  4-  (tip  +  u)  - - -  -1 - - -  =  - 

rit  '  (h  Pp  +  p  (hi  pe  +  p 

fiT  R(T„  +T)  fi( ue  +  u)  fi(Te  +T)  U2  -  (ue  +  u)U, 

__  _)_  +  ( Up  +  u)  — 

lit  cv  fix  fix  (pe  +  p)cv 

In  order  to  proceed  with  easier  notation,  the  perturbation  states  are  denoted  without. 

sign  and  they  should  not  be  mistaken  with  their  state  functions.  Moreover,  the 
nonlinear  terms  are  of  negligible  order  in  the  sulliciently  small  neighborhood  of 

equilibrium  states.  As  an  example,  one  observes 


Pn 


fin 

fix 


fin 

P"e  4—  «  P,,Up 

fix 


fin 

fix 


similar  statements  are  true  about  other  nonlinear  terms,  hence  combining  with  steady 
state  equations  (6-1)  and  (6-3)  one  would  get. 


Pe  dx 


—  p  + 

Pe  dx 


+  RT^  <K  + 

c„  (/X  dx  (he 


(6-22) 


This  set  of  perturbed  equations  are  linear  and  their  characteristics  depend  on  both 
equilibrium  states  and  their  derivatives. 
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7.  LYAPUNOV  STABILITY  FOR  HYPERBOLIC  DISTRIBUTED  PARAMETER  SYSTEMS 


In  this  section  the  Lyapunov  stability  theorem  is  modified  and  applied  to  the  system  of 
partial  differential  equations  of  (6-20)  to  (6-22).  These  equations  represent  the  dynamics  of 
perturbation  state  v  about  an  equilibrium  state  vector  v„.  Hence  by  Lyapunov  method  one  can 
find  whether  anv  deviation  from  the  equilibrium  state  is  stable  or  it  grows  unboundedly  with 
time.  To  be  consistent  with  the  notation  of  sections  3  and  4  one  could  rearrange  equation  (6- 
20)  to  (6-22)  to  the  canonical  form  of  the  evolution  equation.  Here  the  notation  (')  is  used  to 
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(7-1) 

It  is  clear  that  v  -----  0  (i.e..  null  solution)  represents  the  condition  where  there  is  no  variation 
from  the  equilibrium  state. 

The  characteristic  directions  (eigenvalues)  of  this  system  of  equations  are: 


A,  -  u„  ,  A o  ti„  V'-.RTp  ■  A 3  u„  Y/SR Tn  - 

Therefore,  (7-1)  represents  a  system  of  hyperbolic  partial  differential  equations.  For  a  general 
hyperbolic  system  of  the  form 


'v 

-I 


A ' '  v 


B  v  L  v 


with  v  being  a  vector  valued  function  of  dimension  n.  The  following  theorem  indicates  that  the 


operator  L  is  the  generator  of  a  semigroup  on  D|L)  L"  |x<,0./  ,  E"'  . 


theorem  <-l:  1  he  linear  operator  !.  defined  In  (7-.1)  with  A  being  symmetric  is  the  generator 


of  semigroups  (groups  if  the  number  of  positive  eigenvalues  is  equal  to  the  number  of  the 
negative  eigenvalues,  i.e.,  n,  2 )  for  the  solution  of  systems  of  equations  (7-2). 

ft  is  possible  to  show  that  the  system  of  equation  (7-1)  can  be  made  into  the  form  of  (7-2) 
with  a  symmetric  A  matrix. 

For  simplicity  — and  — are  denoted  bv  subscripts  ( •  )t  and  (-)x  respectively.  As  shown 
't  •  'x 

below,  by  dividing  equation  (6-20)  by  (/>„),  (6-21)  by  (A/ RTe  )  and  (6-22)  by  ( N/q  — 1  T,)  the 
following  symmetric  A  matrix  can  be  obtained. 
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it  can  be  shown  that 
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(7-9) 


Construction  of  Lyapunov  functional 


To  construct  a  Lyapunov  functional  it  is  sufficient  to  select  the  functional  as  an  equivalent 
norm  of  a  Hilbert  space  of  states  of  equations  (7-7),  i.c.,  from  the  discussion  in  sections  2  and  4 
this  functional  would  be  a  bilinear  form  as  follows 


=  '  v,  v  '2  =  •  ‘  v,  S  v  —  f  vT  S(x)  v  dx 

o 


(7-10) 


where  S  is  symmetric  positive  definite  and  bounded  linear  operator.  In  order  to  satisfy  the 
sufficiency  condition  of  Zubov’s  theorem  for  the  stability  of  system  of  (7-2)  one  has  to  show 

that  there  exists  a  S(x)  for  which  ~ —  0. 


dt 


Equation  (7-10)  yields: 


I  1  /■  m 

— - —  =  f  v  S(x)  v  dx  4 -  /  vT  S(x)  v  dx 

t  l\  A 


Conjugate  operators  are  defined  based  on  bilinear  form  operation  as  following 


z.Sv  ■  S  z,y  S  —  conjugate  of  S. 


Here,  since  S  is  a  real  operator  S  -  S1 


J  vr  S(x)  v  dx  =  v.Sv  ■  —  ■  S*  v,  v  >  —  j  (STv)T  v  dx 


Therefore, 


— —  -  J  \-T  S(x)  v  dx  +-  J  (ST(x)  v)T  v  dx 
<*t  o  o 


—  ■  v ,  S v  ■  *-  ■  S  rv,  v  > 


The  operator  S(x)  is  chosen  to  be  symmetric  S  =  ST,  then 


=  2  - Sv,  v  •  =  2  f  vT  S(x)  v  dx 


(7-11) 


Now,  from  equations  (7-7),  one  can  substitute  for  v  (i.e.,  vt)  into  (7-11).  This  results  in  the 
following. 


■  S v , (  -  Avx)  -  +  2  -  Sv,(  -  Bv)  • 


I  vT  S(x)  [A  vxj  dx  -2  f  vT  S ( x ) f  13  v|  dx 


where. 
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Using  the  same  discussion  about  conjugate  of  operator  S A  the  following  can  be  derived. 


f  v  1  S(x)  A  — — dx 
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(7-13) 
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(7-14) 


If  [S(x)  Aj  =  S(x)  AjT  equation  (7-13)  becomes  identical  to  (7-14)  and  hence  by  substitution  into 
(7-12),  one  finds 


(  /  (  ^  ( 
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Therefore, 
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2  S(x)B  v  dx 


(7-15) 


Based  on  the  resulting  equation  (7-15),  the  following  conditions  must  be  satisfied  for 
d  J 

guaranteeing  — ; —  ^  t) 

dt 

a)  Boundary  constraint: 

vT(U  )  S(f  )  A(/  )  v(t/  )  -t  vT(t ,0)  S(0)  A ( 0 )  v ( t , 0 )  0  (7-16) 

b)  Interior  constraint: 

f — —  A  •  S(x)  — —  2s(x)  b!  must  be  negative  definite  for  x  (  (0,/^ )  (7-17) 

.  'x  ■  'x  I 


The  operator  A  in  perturbed  linear  model  of  MPD  as  derived  in  equations  (7-6)  to  (7-6)  has 
already  been  put  into  symmetric  form.  Therefore  it  is  sufficient  to  find  a  symmetric  matrix 
function  S(x)  with  discussed  properties  such  that  (7-16)  and  (7-17)  are  satisfied. 

It  should  be  mentioned  that  the  perturbations  from  equilibrium  states  which  are  governed 
by  (7-6)  to  (7-6)  are  subject  to  the  following  initial  and  boundary  conditions 


v|l).\)  V () ( X )  .  v ( t , 0 )  0 


(7-18) 


where 


v0(x) «  L2(oy  ) 


These  conditions  represent  a  situation  where  t He  plasma  Hows  into  the  thruster  with  aflixed 
state  (i.e.,  invariant  with  time).  Due  to  the  perturbation  in  the  control  inputs,  there  is  an 
initial  perturbation  inside  the  thruster  which  is  represented  by  V0(x)  in  (7-18).  Considering 
S(x)  an  operator  even  simpler  than  what,  is  required,  i.e., 


q<*)  0  0 

S(x)  -  0  q(x)  0  =  q(x)  I 

0  0  I 


where  qjx)  is  a  scaler  function,  then  for  S(x)  to  be  positive  definite  it  is  sufficient  to  have, 
qjx)  '  0  for  x  <  10, /"  .  Therefore,  the  inequality  in  condition  (7-10)  would  become 


q( f  )  VT(t,C  )  A(C  )  v(t/)  4-0  0 


where  v(t,0)  -  0. 


It  can  be  seen  that  A(t  )  must  be  positive  definite. 


A(/  )  -  A ( x )  |  ,  ,  VhtT 


l)KT„ 


The  eigenvalues  of  A  are: 


X,  --in  ■  X-.  V'.KT,  ,  X3  !i„  VoHtT 


At  x  -  t  for  accelerating  (low  it  is  practically  desirable  to  have  M  •  1  or  il  vSRT. 


Hence  in  this  case  A(f  )  will  he  positive  definite  and  condition  (7-16)  is  satisfied.  To  evaluate 

,  tS  ■  'A 

the  requirement  for  condition  (7-17).  one  has  to  tind  and  — 


S'  m'IM 

.  'x 


—  A'  -  -VWf.  ~ 


'/-DRT,  — 


-Vl-i-DKT. 

‘  p 


S'fx)  A( x )  *  S|x)  A'(x )  2S( x )  B(x)  q(x)  |  A<x)  A'(x)  -  2  B(x) 

(  q(x) 


Since  (](x)  ■  l)  one  should  have  the  following  condition: 

^  ==Q<*> 

q(x) 

Mi,  =  [Q( x)  A ( x )  +  A'(x) 

where  M[  must  he  negative  definite. 
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where  ax  =  V  R1  „ . 


One  way  to  find  an  apriori  requirement  on  the  system  such  that  stability  of  perturbation  is 
guaranteed,  i.e.,  condition  (7-18)  be  satisfied,  is  to  let  Q(x)  be  negative  and 


|Q(x)|  ’>  max  f 


f-  |3/2~l  -I— I  +  2  1^1^2  1-^1 

i„  up  R 1  e 


hence 


ueQ  ilxQ  0 

Ml  =  axQ  ueQ  V~/— 1  axQ 

0  V-,-1  aTQ  ueQ 


Then  condition  (7-18)  would  become 


u„  aT 


-  aT  un  V~  -1  ax  is  positive  definite 


Again. 


X,  =  u„  ,  X2  ----  t-  VoRTp  .  X3  -  Up  \ArTp 


for  -  to  be  positive  definite,  one  should  have  X3  •  0 

Q 


u„  •  “\/  I R Tp  ,  M  I  for  x  (  (0,1) 

This  restates  the  situation  that.  Ilow  is  supersonic  through  the  interior  of  the  MPD  thruster. 
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conclude  tlvat  a  su 


„kiclll  ^  f*  *.  "  lh' ayst"" 


Th""0re  '  .  (h  ,  lh,  Uo.  ,eSi,„e  remains  supersonic  tbcou6bnut  tbe 

respect  to  initial  perturbations  «  that 


thruster. 
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The  analysis  of  Section  7  indicates  that  in  the  case  of  supersonic  thrusters  there 
was  sufficient  conditions  for  system  stability  after  an  initial  perturbation.  However,  in 
general,  one  would  like  to  have  the  system  stabilized  for  all  How  patterns  (all  steady 
state  solutions),  which  are  of  interest  as  discussed  in  Section  6.  In  order  to  address 
this  problem  the  inherent  characteristics  of  the  system  should  be  analyzed.  Such 
analysis  often  involves  information  about  characteristics  of  the  system  of  partial 
differential  equations. 

Recalling  the  system  of  equations  (7-7)  for  a  general  case 

+  A  ~~  +  Bv  =  0.  (8-1) 

(H  <7X 

Where  A(x)  is  symmetric  and  has  eigenvalues  of  the  general  form; 

X,(x)  <  \2(x)  < .  <  Xp(x)  <  0  <  Xp+1(x)  <  •  •  •  <  Xn(x). 

The  eigenvalues  show  the  directions  of  the  characteristic  lines  of  PDE.  Also  there 
exists  a  continuously  differentiable  matrix  O(x),  based  on  the  system  eigenvectors  such 
that  0-1(x)  A(x)  O(x)  =  \(x).  If  one  considers  a  new  set  of  states  w  such  that 


v  =  0(x)w 

then  substitution  of  v  into  (8-1)  results  in 


O(x)^-  +  A(x)Q(x)—  + 


A(x)^^-  -j-  BO(x) 
(h. 


w  =  0 


or 


^  +  0_1(x)A(x)0(x)-~  +  O1 

Ht  <h c 


A(x)  +  DO(x) 

(h. 


w  =  0 


Therefore  one  can  conclude; 


Since  A(x)  is  a  diagonal  matrix  of  the  system  eigenvalues  and  from  the  fact  that  X,  to 
Xp  are  negative  and  Xp+1  to  Xn  are  positive,  then  A  can  be  decomposed  as  follows. 


A  = 


A  =  diagfX,,  •  •  •  ,Xp) 

A  +  =  diag(Xp  +  I . Xn) 

The  corresponding  decomposition  can  be  applied  to  the  states 


A  O 
0  A * 


(8-3) 


w  = 


(8-1) 


Therefore  (8-2)  can  be  reduced  to  a  set  of  ordinary  differential  equations  of  the 
following  form. 


dt 


_ 

ih  +  k  ihi 


<'Kvk  <9wk 


+ 


dx 

dt 


<)wk 

ih. 


=  -  /4(x)w 


ih 

k  =  1,  2.  •  •  ‘  ,n  (8-5) 

where  wk  is  kth  element  of  the  vector  valued  function  w.  Moreover,  Xk  represents  the 

direction  of  the  kth  characteristic  line,  hence  ^ d  is  the  directional  derivative  along 

dt 

the  corresponding  characteristic  line.  On  the  right  hand  side  of  (8-5)  fik,  is  the  kth 
row  of  the  matrix  valued  function  /f(x).  The  set  of  n  ordinary  differential  equations 
(8-5),  are  coupled  by  the  term  /7k(x)w.  Initial  values  of  t  0  can  be  given  as 


w(x.O)  =  wG  <(L21o/  ]:  K“) 

Considering  i  he  boundary  conditions,  at  a  point  on  the  boundary  x  0.  t  -i  as 
shown  in  figure  (8-1),  one  finds  that  the  characteristics  with  negative  and  positive 
eigenvalues  arrive  at  that  point  with  negative  and  positive  slopes  respectively.  The 
incoming  information  consists  of  values  of  wk  associated  with  characterist ic  line 

Ck(0.  tn)  with  negative  slope  for  k  1 .  p.  The  "outgoing  information"  consists  of 

values  of  wk  associated  with  positive  slope  characteristic  Ck(o,  tn),  k  p*  I .  n. 

Hence  along  the  boundary  x  0  the  values  of  w‘  should  be  known.  It  is  clear  that 
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along  boundary  x  =  1,  the  orientation  of  characteristics  will  be  reversed  and  hence 
the  values  of  w~  should  be  known. 


At  this  stage  one  can  impose  the  question  of  stabilizability  of  linear  hyperbolic 
systems,  namely  those  systems  such  as  MPD  thruster  with  an  initial  perturbation. 
The  stabilizability  question  would  be  as  such:  what  control  inputs  should  be  used  to 
bring  the  system  to  a  stable  state?  This  question  can  be  approached  in  two  different 
ways.  The  first  approach  is  the  idea  of  null  controllability  which  is  more  restrictive 
and  considers  determination  of  a  control  (applied  on  the  boundary)  such  that  for  a 
given  initial  perturbation  w0,  after  time  t  >  0, 


w(x.T)  =  0. 

It  can  be  seen  that  the  required  time  T  should  satisfy  the  following  inequality 


T  >-/  -Ji-  +  / 


dx 


Xp(x)  o  Vl(X) 

By  this  approach  t he  time-space  plane  (domain  of  the  system)  can  be  divided  into 
three  separate  regions  as  shown  in  figure  (8-2).  After  finding  the  values  of  states  in 
each  region,  the  value  of  w(f  ,t)  will  be  evaluated.  Since 


w(f  ,t)  = 


w~(f  ,t) 

w  +  (F  ,t ) 


then  knowledge  of  both  incoming  and  outgoing  information  along  x  =  /  would  lead  to 


the  determination  of  a  boundary  value  control  at  this  boundary.  20 


1  he  second  approach  is  more  or  less  based  on  the  similar  ideas  as  those  of  Section 
7.  In  this  approach  an  energy  function  (i.e.  Lyapunov  function)  is  defined.  Depending 
on  the  type  of  system,  two  types  of  controllers  can  be  sought.  These  are  boundary 
controller  or  body  force  controller.  For  a  wave  equation  (linear  and  nonlinear)  it  is 
shown  la  that  a  control  force  at  the  boundary  ran  steer  the  system  to  a  stablized 


V- .V.'.-.V.V’v  V.v.v..-.  -Y--V.V ,  .  .  «Y  .■ 


state.  Whereas  in  the  case  of  a  linear  wave  equation  with  a  "body  force"  (i.e.  internal 
controller)  it  is  shown  { 1 6 j  that  the  system  can  be  controlled  to  a  "zero  energy"  state. 
In  both  examples  the  energy  functions  were  proposed  as: 

E  =  —J  ft  ’  +  C’  ft  ’dx 
2  \  (h  dx 

Although  both  of  these  controllability  approaches  have  been  applied  to  much 
simpler  systems,  it  is  desirable  to  answer  some  of  these  controllability  questions  in 
reference  to  the  MPD  system  via  these  approaches  in  the  future  research. 
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■ABSTRACT 

A  magneto-plasma-dynamic  thruster  is  subject  to  a  number  of  instabilities,  some  of 
them  associated  with  ionization  or  electro-thermal  instability.  An  attempt  has  been 
made  (Part  I)  to  establish  the  conditions  under  which  imbalances  in  power  input  and 
electron  diffusion  can  lead  to  instability.  In  Part  II  a  new  approach  to  stability  of 
MPD  engines  is  presented.  This  approach  is  based  on  extension  of  the  Lyapunov 
direct  stability  theorem  to  distributed  parameter  systems.  Problem  formulation  in  an 
operator  form  is  presented.  The  Lyapunov  function  is  extended  to  a  Lyapunov 
functional  which  is  an  integration  over  space.  This  approach  has  been  applied  to  a 
magneto-gas  dynamic  problem,  namely  a  simplified  form  of  an  MPD  engine.  The 
stability  results  are  presented  and  discussed.  Wave  characteristics 'in  both  transverse 
and  longitudinal  modes  are  discussed,  and  upper  bounds  for  wave  speeds  based  on 
Alfven  speed  are  derived. 
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Background 

A  MPD  thruster,  shown  schematically  in  Figure  1,  may  be  divided  into  three 
regions:  (l)  initial  ionization  region  with  heating  of  electrons,  rapid  dissociation  and 
ionization  of  heavy  particles  and  electron  pressure  large  compared  to  the  partial 
pressure  of  heavy  particles;  (2)  uniform  arc  region  with  electron  temperature  nearly  in 
equilibrium,  heavy  particles  being  heated  by  transfer  of  energy  from  electrons 
undergoing  Joule  heating  and  electron  temperature  and  current  density  nearly 
constant;  and  (3)  nonuniform  arc  region  with  electron  properties  varying  appreciably 
while  the  advection  velocity  is  transonic  or  greater. 

One  important  factor  in  region  (2)  is  that  the  Joule  heating  undergoes  a  change 
when  the  plasma  current  attains  a  constant  value  in  space.  Thus,  initially,  Joule 
heating  can  be  expressed  by  the  relation,  namely 

JH(Te)  a  Ip2  T^2  (1) 

while,  when  the  current  attains  the  constant  value,  it  becomes  the  following. 

JH(Te)  a  V2T3/2  (2) 

In  Eqs.  (1)  and  (2),  Ip,Te  and  V  represent  the  plasma  current,  electron  temperature 
and  voltage,  respectively.  It  is  important  to  note  that  a  characteristic  diffusion  time 
has  to  elapse  before  Joule  heating  changes  from  that  given  by  Eq.  (1)  to  that  by  Eq. 
(2).  At  the  end  of  such  diffusion  time,  one  can  equate  JH  in  Eqs.  (1)  and  (2)  and 
obtain  an  expression  for  the  voltage,  namely 

V  a  Ip/Te3/2  .  (3) 

Our  main  objective  is  to  examine  the  structure  of  plasma,  as  given  by  "  set  of 
describing  equations,  and  to  determine  if  there  arises  a  branching  and  the  consequence 
when  such  a  state  of  constant  current  arises. 


Describing  Equations 


The  partially  ionized  gas  is  considered  u  ider  the  hydrodynamic  approximation, 
given  by  Refs.  9-10.  While  the  resulting  equations  are  nonlinear  partial  differential 
equations,  we  consider  the  lowest  order  approximation  to  those  in  terms  of  the 
following  ordinary  differential  equations.  They  will  be  adequate  for  illustrating  the 
effects  of  the  plasma  current  attaining  a  constant  value  over  a  characteristic  period  of 
time.  The  equations  of  interest  are  for  the  number  density,  temperature  and  energy 
balance  of  relevant  species. 


dNe 

dt 


a 


(4) 


From  Eq.  (9),  one  can  solve  for  T,  and  substitute  for  the  Joule  heating  term  from  Eq. 
(1).  Similarly  Eq.  (10)  provides  another  relation  between  T,  and  Te.  It  is  then  of 
interest  to  examine  the  two  curves  when  Te  is  small  and  when  it  is  sufficiently  large. 
After  some  algebra,  the  two  curves  can  be  utilized  to  deduce  (Ref.  11)  the  domains  in 
which  (dTe/dt)  and  (dT,/dt)  are  positive  or  negative,  as  shown  in  Fig.  2.  The  point  of 
intersection  of  the  curves  is  a  stable  equilibrium  point. 

Referring  to  Eqs.  (9)  and  (10),  it  is  next  necessary  to  substitute  for  Joule  heating 
from  Eq.  (2).  For  Te  — *  0,4>e  starts  at  the  origin.  When  Te  becomes  finite,  the  relation 
between  T;  and  Te  depends  upon  the  diffusion  coefficients  SDe  and  S^j.  In  particular, 
several  solutions  can  be  obtained  depending  upon  the  enhancement  of  diffusion  in  the 
gradient  zone,  that  is,  depending  upon  the  value  of  a.  Thus,  while  for  small  Te  the 
value  of  T,  is  nearly  equal  to  Te,  at  larger  values  of  Te,Tj  decreases  and  becomes 
negative,  the  change  depending  upon  the  extent  of  enhanced  diffusion  or  the  value  of 
a.  Several  illustrative  cases  can  be  formulated. 

(a)  Case  of  small  enhancement  of  diffusion:  Two  subcases  can  be  distinguished 
when  ot  is  small,  that  is,  less  than  a  certain  value  that  is  a  function  of  V,Te  and  Ne. 
The  first  subcase  is  obtained  for  small  input  power.  Then,  as  shown  in  Fig.  3,  the 
system  returns  to  the  origin  or,  for  large  Te,  a  runaway  type  structure  is  obtained. 
The  point  of  intersection  of  the  (dTe/dt)  and  (dT,/dt)  curves  is  then  a  saddle  point. 

The  second  subcase  pertains  to  large  input  power.  Then,  there  may  arise, 
depending  upon  the  values  of  system  parameters,  one  of  two  possibilities:  (1)  As  in 
Fig.  4,  the  origin  may  become  unstable  while  the  two  trajectories  of  the  solutions  do 
not  interest  or  (2)  as  in  Fig.  5,  there  may  arise  two  points  of  intersection.  In  the  latter 
case,  both  the  origin  and  the  second  point  of  intersection  are  unstable  while  the  first 
point  of  intersection  is  a  stable  node. 

(b)  Case  of  large  enhancement  of  diffusion:  In  the  case  at  is  large,  the  parameters 
of  influence  are  Ne  and  the  power  input.  For  a  given  value  of  power  input,  there  is  a 
value  of  Ne  below  which  a  stable  node  is  obtained  as  shown  in  Fig.  6.  Above  such  a 
value  of  Ne,  the  two  trajectories  do  not  interest,  as  shown  in  Fig  7.  In  other  words,  at 
any  value  of  Ne  there  is  a  maximum  amount  of  power  that  can  be  applied,  beyond 
which  value  of  power  the  enhanced  diffusion  will  drive  the  system  towards  the  origin. 

Disc  ussio  n 

In  the  region  where  plasma  is  being  supplied  with  energy  for  the  heating  of  ions 
through  the  Joule  heating  of  electrons,  a  branching  in  the  structure  of  plasma  has 
been  shown  to  arise  when  the  current  across  the  plasma  becomes  constant  and  remains 
at  that  value  over  a  characteristic  period  of  time  related  to  diffusion  of  electrons.  The 
diffusion  becomes  enhanced  in  the  gradient  zone.  When  the  net  diffusion  is  small,  the 
various  possible  states  of  plasma  depend  upon  the  amount  of  power  supplied.  On  the 
other  hand,  when  the  net  diffusion  is  large,  there  exists  a  number  density  of  electrons 
at  each  given  value  of  power  input  above  which  there  appears  a  disruption  in  the 
structure  of  plasma.  This  provides  a  basis  for  regulating  power  input  in  the 
equilibrium  region.  The  input  parameters  for  such  a  control  system  become  the 
electron  number  density  and  ion  temperature. 


Problem  Formulation; 


Consider  the  MPD  engine  shown  in  Figure  1.  The  plasma  dynamic  equations  for 
this  engine  consist  of  Maxwell’s  equations,  Ohm’s  law,  conservation  of  electric  charge, 
equation  of  state  (ideal  gas  law)  and  a  set  of  mass,  momentum  and  energy  equations 
[24).  These  equations  can  be  written  as  follows: 

Maxwell  equations: 
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Equation  of  state  (ideal  gas  law): 


P  =  RpT 


Conservation  of  Mass: 
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current  density  J  and  net  electric  charge  pe  are 
J  =  J  (x,t),  pe  =  pe(x, t)  . 


The  one-dimensional  assumption  results  in:  =  0,  =  0, 
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Based  on  the  aforementioned  treatment  of  the  problem,  the  following  describing 
equations  can  be  derived. 


Maxwell’s  equations: 
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Generalized  Ohm's  law: 
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Conservation  of  electric  charge: 
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component  of  speed  of  Alfven  wave. 


(ii.)  Longitudinal  Mode 

The  state  equation  for  this  mode  can  be  reduced  to 
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where  p  =  -c— ,  T  =  -^7-,  Vv  =  \  /  —  Hy.  The  parameter  Vv  is  defined  as 


Po  lo  '  V  Po 

y-component  of  the  speed  of  Alfven  wave, 

Lyapunov  Functional  and  Stability  Analysis 

In  this  section  the  Lyapunov  Functional  approach  is  applied  to  each  mode  of  the 
plasma  dynamics.  The  stability  results  are  derived  and  discussed. 

(i.)  For  the  transverse  mode  the  general  operator  form  of  the  equation  (17)  and  (18) 
evolution  equation  would  be 


In  order  to  have  T  and  X  functions  independent  from  x  and  t,  respectively,  it  is 

T  X"  T  X  X" 

required  that  —  =*  f(x,t),  *  f(x,t).  From  — =  const,  and^—  =  const 
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it  is  conclusive  to  represent  X  functions  in  terms  of  a  real  periodic  function 
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where  A  =  complex  conjugate  of  A 
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The  point  spectrum  of  operator  A,  i.e.,  <rp(A)  can  be  found 
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the  following  can  be  resulted. 
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In  this  case  (28)  shows  that  V  is  negative  definite  and  that  proves  the  stability  of 
the  transverse  mode  using  Lyapunov  approach. 


(ii.)  For  longitudinal  mode  of  wave  propagation  the  generic  form  of  evolution  equation 
(24)  is  considered 
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and  Z(o,t)  =  Z(f  ,t)  —  0  . 


In  this  case  A  in  equation  (24)  is  a  linear  operator  with  the  domain  in  a  separable 
Hilbert  space  of  state  function  which  maps  Z  onto  itself.  From  equations  (19)  to 
(23)  A  is  formed  as, 
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In  general  Sn  can  be  written  as  a  combination  of  real  and  imaginary  parts. 


Sa  =  Re(Sn)  +  i  Im(Sn)  . 


This  resulted  in 
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Obviously,  in  order  to  have  a  wave  (under  damped  condition)  u>2  has  to  be  positive. 


Otherwise,  (o>2  <  0)  there  is  no  wave  propagation  in  the  transverse  mode.  Hence, 


wave  speed  can  be  defined  as 
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which  has  an  elliptical  shape,  as  shown  in  Figure  8.  The  maximum  of  the  transverse 
wave  speed  is  limited  by  Alfven  speed  in  the  x-direction  in  the  case  of  iqj  =  u. 
Moveover,  for  the  transverse  wave  to  exist,  the  maximum  difference  in  —  u) 
caused  by  minimum  Xn  (i.e.  Xn  rain  =  n).  Therefore,  if 
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the  generation  of  transverse  wave  is  plausible. 


For  the  longitudinal  waves,  due  to  the  complexity  of  characteristic  equation  the 
case  that  u  —  K  =  0  and  #  0  is  considered.  The  characteristic  equation,  then  will 
be, 
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here  V2  =  V2  +  V2.  For  the  case  of  approaching  infinity  S3  is  zero  and 
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u,  v,  w:  velocity  vector  components  in  x,y  and  z  directions 

x,  y,  z:  Spatial  coordinates 

X.:  Spatial  coordinates  for  j  =  1,  2,  3 


Greeks 


Dielectric  constant 


/i,  v  Viscosity,  kinematic  viscosity 

He:  Magnetic  permeability 

i 

p,  p:  Plasma  density,  perturbation  in  density 

pe:  Excess  electric  charge/volume 

cr.  Electrical  conductivity 


subscripts 

x,  y,  z:  Variable  subscribed  in  the  direction  of  x,  y  and  z. 

Notation 

Tjj,  Tw:  Time  dependent  part  of  variable  in  subscript 
Xh,  Xw:  x  dependent  part  of  variable  in  subscript. 
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Abstract 


An  attempt  has  been  made  to  set  up  a  model  for  heat  transfer  to  the  cathode  of  a 
MPD  thruster  and  to  conduct  an  observability  analysis  of  a  model  energy  balance 
equation.  Considering  a  monoatomic  gas  that  is  singly  ionized  and  is  in  a  state  of 
nonequilibrium  with  different  species  at  different  temperatures,  the  energy  balance  in 
the  vicinity  of  the  electrode  is  described  in  terms  of  seven  regions  including  a  "mushy" 
region  involving  phase  change  processes.  The  observability  analysis  is  conducted  on  a 
model  energy  balance  equation  applicable  to  one  of  the  regions.  Since  a  MPD  engine  is 
a  nonlinear  distributed  parameter  system,  the  concept  of  observability  analysis  for 
such  systems  is  derived  and  presented.  An  observability  criteria  for  finite  and  infinite 
dimensional  systems,  based  on  the  transformation  of  nonlinear  systems  to  observable 
canonical  form,  is  studied.  A  less  restrictive  approach  for  observability  of  distributed 
parameter  systems  is  presented.  This  technique  is  applied  to  two  cases  of  available 
measurements  of  the  MPD  engine.  The  resulting  conditions  on  the  control  inputs  to 
the  engine  are  presented. 


This  research  has  been  supported  under  AFOSR  Grants  85-0335  and  80-0278,  with  Dr.  J.  Tishkoff  as 
Technical  Monitor  of  the  project. 


INTRODUCTION 


The  electrode  region  of  a  magneto-plasma-dynamic  (MPD)  engine  continues  to 
present  several  complexities  both  in  analysis  and  measurements  that  can  serve  to 
explain  phenomenology  and  to  permit  predictions  of  such  quantities  as  heat  transfer 
and  wear  [1-6].  The  complexities  are  further  accentuated  when  the  operating 
conditions  are  close  to  the  occurrence  of  onset  instability  and  other  critical  regimes  [7- 
9].  The  electrode  processes  are  of  interest  both  during  the  start-up  phase  of 
continuously  operated  systems  and  in  pulsed  devices.  In  all  cases  attention  has  to  be 
paid  to  the  state  of  the  plasma  and  the  processes  that  are  dominant  in  different  zones 
as  the  electrode  surface  is  approached  from  the  free  stream  condition.  In  general,  the 
plasma  may  be  multi-component,  multi-temperature,  including  elastic  and  inelastic 
collisions,  with  complex  reaction  rates,  subjected  to  Joule  heating,  and  acted  on  by 
induced  and  applied  magnetic  fields.  Furthermore,  one  has  to  account  for  electrode 
surface  reactions  and  resulting  changes  both  in  the  plasma  and  the  mechanical 
structure  of  the  electrode  itself.  A  number  of  zones  can  be  identified  in  the  vicinity  of 
the  electrode  based  on  characteristic  length  scales  such  as  the  recombination  length 
(^),  the  molecular  mean  free  path  (4,e,i,  a>  e  and  i  referring  to  atom,  electron  and  ion, 
respectively)  and  the  Debye  length  (<q).  The  rate-governed  processes  in  those  zones 
will  depend  further  upon  characteristic  times. 

A  beginning  has  been  made  in  Reference  [10]  on  the  analysis  of  a  plasma  boundary 
layer  under  the  assumption  of  an  equilibrium  plasma  in  the  absence  of  Joule  heating 
and  any  ambient  magnetic  field.  Although  the  boundary  layer  pertains  to  a  power 
generator  (and  not  a  plasma  generation  device),  important  findings  have  been  made  on 
the  relative  importance  of  inertia  of  electrons  and  ions  and  of  charge  separation  in  the 
vicinity  of  electrodes  under  various  levels  of  current  input.  The  energy  balance 
equation  is  obviously  left  out  of  account  in  the  analysis.  On  the  other  hand,  inclusion 
of  considerations  of  nonequilibrium  and  Joule  heating  make  it  imperative  to  include 
the  energy  balance  equation  in  the  analysis. 

The  analysis  of  a  complex  flowfield  such  as  that  of  a  MPD  engine  can  be 
undertaken  from  several  points  of  view:  prediction  of  performance,  analysis  for 
distinguishing  and  ordering  the  importance  of  various  processes  and  determination  of 
stability,  optimality,  observability  and  controllability  of  the  system  represented  by  a 
set  of  describing  equations.  The  latter  has  been  the  basis  of  recent  investigations 
undertaken  at  Purdue  University. 

The  objectives  of  the  current  work  in  that  context  are  two:  (1)  To  obtain  a 
reasonably  detailed  set  of  describing  equations  for  the  vicinity  of  the  electrode, 
including  Joule  heating  and  energy  balance,  and  (2)  to  present  the  concept  of 
observability  for  such  nonlinear  systems  as  MPD  engines,  where  system  parameters 
change  with  both  time  and  space  (Distributed  Parameter  System  -  DPS). 

In  general,  distributed  parameter  systems  are  represented  by  partial  differential 
equations,  relating  informations  on  spatial  domain  of  the  system  to  the  time  domain. 
These  systems  can  be  considered  infinite  dimensional  as  opposed  to  finite  dimensional 
notion  for  lumped  parameter  systems. 

Control  of  the  DPS's  are  usually  achieved  by  a  parameter  approximation  of  the 
system.  This  approach,  using  finite  element  or  finite  difference  technique  has  been 


useful  for  control  design  of  some  systems.  In  general,  intrinsic  properties  of  DPS  can 
not  be  predicted.  The  stability,  controllability  and  observability  considerations  of 
such  systems  often  can  not  be  achieved  with  discretized  models.  The  formulation  and 
treatment  of  continuous  model  approach  from  the  standpoint  of  control  theory  has 
been  investigated  [11,12]. 

Many  researchers  have  studied  the  stability  of  linear  DPS  [13,14].  One  of  the  first 
attempts  using  Lyapunov’s  direct  method  for  dimensional  systems  was  made  by  [15]. 
Based  on  an  abstract  theory  of  Lyapunov’s  stability  scheme,  references  [18,17]  have 
defined  Lyapunov’s  functions  on  Hilbert  and  Banach  spaces.  An  application  of  this 
approach  has  been  studied  for  magneto-plasma-dynamic  (MPD)  system  modeled  as  a 
DPS  [18]. 

The  controllability  of  linear  DPS  has  been  formulated  and  studied  by  [12,14,19]. 
Also  the  researchers  [20]  have  investigated  and  derived  a  controllability  criteria  for 
special  classes  of  DPS.  Reference  [21]  is  attempting  to  relate  the  controllability  of 
nonlinear  finite  dimensional  systems  originated  by  [22,23]  to  a  class  of  nonlinear  DPS. 

In  a  distributed  parameter  system  the  question  of  observability  may  cover  two 
types  of  information:  the  criteria  which  determines  whether  the  measurements  contain 
sufficient  conditions  to  uniquely  describe  the  state  of  the  system;  and  the  location  of 
sensors  in  order  to  provide  feedback  information.  This  paper  is  primarily  devoted  to 
the  first  aspect  of  the  problem  for  the  purpose  of  establishing  a  criteria  for  the 
observability  of  DPS.  One  of  the  initial  investigations  in  this  area,  Wang  [22],  has 
discussed  certain  aspects  of  observability  for  distributed  parameter  systems  and  has 
given  a  definition  with  respect  to  initial  state  recovery  of  systems.  Goodson  and  Klein 
[24]  define  observability  as  the  ability  to  establish  the  uniqueness  of  a  solution  of  the 
system.  Other  investigators  [25-29]  have  studied  the  observability  of  systems  which 
belong  to  the  classes  of  linear  distributed  parameter  systems.  For  such  systems,  often, 
it  is  necessary  and  sufficient  to  have  a  non-zero  inner  product  of  measurement 
distribution  with  infinite  dimensional  eigenvector  of  the  system’s  differential  operator. 

For  nonlinear  systems,  in  general,  the  eigenfunction  can  not  be  determined  and 
consequently  the  methods  used  for  linear  systems  would  not  be  applicable  to  nonlinear 
infinite  dimensional  systems.  For  nonlinear  finite  dimensional  systems  many 
researchers  have  tried  to  transform  the  nonlinear  system  equations  into  a  set  of  linear 
equations  with  a  specific  canonical  form.  References  [30,31]  present  the  conditions  on 
the  existence  of  such  transformations.  An  observable  canonical  form  for  nonlinear 
systems  is  constructed  in  [32]. 

In  this  study,  the  objective  is  to  find  a  region  in  the  state  space  for  which  the 
measurements  of  the  states  would  lead  to  reconstruction  of  the  system  ( i . e . ,  the  unique 
behavior  of  the  system  is  feasible). 

GENERIC  MODEL  FOR  ELECTRODE  REGION 

Figure  1  presents  a  schematic  of  the  electrode  region  of  a  MPD  engine.  A  thermal 
boundary  layer  is  postulated  that  may  be  different  from  the  momentum  boundary 
layer.  The  Debye  length  and  the  mean  free  path  of  species  are  both  small  compared 
to  the  recombination  length  but  may  be  either  less  than  or  be  comparable  to  ^ei. 
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The  working  fluid  is  assumed  to  be  a  monatomic  gas  that  is  undergoing  single 
ionization  so  that  the  species  considered  are  atoms,  ions  and  electrons  obeying 
Maxwell-Boltzmann  statistics.  The  plasma  in  the  free  stream  is  assumed  to  be 
partially  ionized  with  a  possibility  of  both  thermal  and  charge  nonequilibrium.  The 
electron  and  ion  temperatures  (TgT,)  are  therefore  assumed  to  be  different.  The  Saha 
equation  for  reaction  is  modified  in  order  to  take  into  account  non-equilibrium  and 
inelastic  collisions.  The  plasma  is  assumed  to  be  subject  to  Joule  heating  and  to 
induced  and  applied  electro-magnetic  fields. 

The  state  of  plasma  within  and  therefore  may  only  be  described  in  terms  of 
collisionless  particles.  It  may  be  pointed  out  that  ^ei  and  are  not  sharply  defined 
boundaries. 

Within  the  electrode,  two  regions  have  been  identified  in  Figure  1:  the  first  is  a 
region  where  melting  and  evaporation  may  be  taking  place  and  the  second,  a  region 
with  pure  conductive  heat  transfer.  The  first  is  referred  to  as  a  "mushy"  region. 

The  describing  equations  applicable  to  the  various  regions  are  provided  in 
Appendix  I.  At  this  stage  it  is  not  the  objective  either  to  match  the  different  regions 
on  an  asymptotic  basis  or  to  undertake  numerical  calculations.  Hence  no  boundary 
conditions  are  given  for  the  different  regions. 

ENERGY  BALANCE  AT  A  CATHODE 

As  a  particular  example  of  electrode-associated  processes,  we  examine  the  energy 
balance  in  the  vicinity  of  a  cathode.  The  equations  describing  the  state  of  the  plasma 
including  non-equilibrium  and  multi-temperature  effects  are  presented  in  Appendix  I. 
Referring  to  Figure  1,  the  energy  balance  equations  applicable  to  different  regions  are 
presented  in  Appendix  II  along  with  the  associated  current  flux  equation.  It  is 
assumed  that  the  cathode  material  is  catalytic  and  recombination  processes  within  the 
material  give  rise  to  a  flux  of  neutrals  into  the  Debye  region  and  also  to  heat 
generation  within  the  electrode.  A  region,  referred  to  as  the  "mushy"  region,  is 
identified  in  the  electrode  material  adjoining  the  plasma.  In  that  region  it  is  assumed 
that  phase  change  processes,  solid  to  liquid,  liquid  to  vapor  and  also  solid  to  vapor,  are 
expected  to  occur. 

The  energy  balance  equation  at  the  internal  surface  of  the  cathode  marked  in 
Figure  I  as  y=0  may  be  written  as  follows. 

^boiid  Qp  ^Irad  Qsur  |Q  "b  qcond  ~b  Select!  mushy 

where  qp  is  the  net  energy  transfer  to  the  surface  from  the  plasma  due  to  the  kinetic 
energy  and  potential  energy  of  the  particles,  qrad,  the  net  radiation  to  the  surface  from 
the  plasma,  Qj^^the  energy  generation  at  the  surface  due  to  recombination  of  ions 
and  electrons,  Q  ,  is  the  energy  associated  with  phase  change  and,  qcond,  energy 
conducted  into  the  mushy  region. 
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observability  analysis 


The  general  form  for  the  3et  of  governing  equations  of  the  system  under 
investigation  can  be  expressed  as 


. » 


(1) 


y=h(V) 

where  V  is  the  distributed  state  vector  and  y  is  the  scalar  output  (measurement). 

A  nonlinear  transformation  T:V  —►V  with, 

V  =  T(V*)  (2) 

such  that  system  of  (1)  is  transformed  into  a  Brunovsky  canonical  form  [22]  given  by 

V  *  =  AV  *  +  7(y,  U)  =  f*(V‘,  U) 

y  =  CV‘ 


with 


A  = 


0 

0  .  . 

.  0 

0 

1 

0  .  . 

0 

0 

I  .  . 

. 

. 

0 

0  .  . 

• 

• 

0 

, 

... 

0 

0 

0 

0  .  . 

I 

0 

and 


C  =  [0,  0,  ....  0,  1] 


Hence 


V  =  £L(V-)  -  <VT(V)  ■  2L  > 


f(V,  U)  =  <  VT(V')  ■  f*(v*'  U)  > 
* 

by  differentiation  with  respect  to  Vk, 


illYJJl  _  +  VT(V-, . 

avk  avk 


L.H.S.  -  M.  msl 

■>v  avk 


First  term  of  the  R.H.S.  i 


is 


(3) 


(4) 


(5) 


(6) 

(7) 


(8) 

(8-a) 
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F 


j2I_ 

cTVj 


•  VT(V**  •  f*  = 


or 

dwi 


dV 


f(V) 


(8-b) 


dt 


According  to  the  chosen  canonical  form,  - r  for  k  =  1,  n  -1  is 

<3Vk 


0 


Of 

3Vl 


Namely,  (k+l)th  element  is  unity.  Hence  the  second  term  on  the  right  hand  side 

dT 


is 


dT  i)( 


for  1  <  k  <  n— 1 


(8-c) 


ay  ,tv. 


for  k  =  n 


Therefore,  from  (8)  using  8-a,  b,  c,  for  1  <  k  <  n— 1,  the  following  is  obtained: 

dT 


a 


'T 


dt  dT 


cN\. 


'  <-1  W  dVy 


av 


f  = 


[f 

*  f  _  $ 

ii 

ad1  f, 

?y k. 

dV, 

(9) 


where 


f, 


a  t 


ay. 


dT 


is  Lie  bracket  of  the  two  vectors,  f  and  - r-  From  (3),  if  the 

avk 


measurement  y  is  differentiated,  then 


which  can  be  reduced  to 


:iY  _  ^  ™ 

<  n'  ay  r)y' 

<  Vh(V) ,  > 

ay 


=  c 


=  c 


(10) 


(11) 


Elements  of  the  left  hand  side  vector  are 


<  Vh(V)  ,  -^r  >  =  <  Vh(V)  ,  i,~  > 

uV 2  l 


From  Leibritz  formula 

<  Vh(V)  ,  f,  >  =  <  V  <  Vh,  f  >,  >  ~  <  V  <  Vh,  >,  1 

OV  2  Uy  |  CTv  ^ 

But  considering  the  zero  inner  product  in  (12) 

,rr  rtT 

<  Vh,  f,  —b  >  =  <  v  <Vh,  f  >  ,  > 

<w,  avj 

where  <  Vh,  f  >  =  Lf(h)  is  Lie  derivative  of  h  with  respect  to  vector  field  f. 


Therefore, 


<  Vh,  >  =  <  VLf(h),  > 
<)v2  OV , 


Similarly,  it  will  be  found  that 


<  Vh,-— >  -  <  VLf2(h),  > 

'  ll  r  I  \  /  r 


<  Vh,  >  =  <  VLfn-I(h),  — > 

'nr  1  '  '  •  nr 


Oh  >fV  _ 

i-N  ;n 


Vh 

VLf(h) 

VLfZ(h) 

VLf3(h) 


VL“-l(h) 


Lie  derivatives  of  higher  orders  L2(h),  Lf3(h),  .  .  .  are  defined  as 
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L,!(h)  -  Lf(Lf{h)) 


L,k(h)  =  MM . L,(h))  ■  •  •  ) 


In  order  to  have  system  of  equations  (1)  transformable  to  observable  form  of 
equation  (3)  the  transformation  set  T  and  its  derivative  must  exist.  Hence,  - r  must 

v  '  dVv 

exist  and  in  equation  (16)  the  observability  matrix  O  should  be  invertible. 


O  = 


This  condition  can  be  met  by 


Vh 

VLf(h) 


VL,-‘(h) 


det  O  #  0 


(18) 


From  this  point  on  two  questions  are  considered.  One  being  the  condition  for  which 
the  set  of  transformations  T  can  be  found,  i.e.,  solvability  of  partial  differential 
equation  in  terms  of  derivatives  of  T  with  respect  tc  transformed  state  variables,  i.e., 
V,  ,  i=l,  ...,  n.  The  condition  for  solvability  of  this  partial  differential  equation  is 
given  by  Forbenius  solvability  condition,  which  requires  the  set  of  vectors 

r)  T 


<JV 


ad°f, 

adlf,  --^r 

J  ••••) 

ada-lf,  — 

.1  Wi 

.  ^V1 

Wi 

be  involutive  [33).  A  set  of  vectors 
any  two  is  a  linear  combination  of  the 

f,(V),  f,(V) 


i=l,  ...,  n  is  involutive  if  the  Lie  bracket  of 


f; 


t,  i-e. 


=  ck  fk(V) 

k  =  l 


(19) 


This  condition  is  the  tailored  form  of  the  fact  that  each  transformed  state  can  be 
solved  in  terms  of  state  variable 


V,  =  Vj’fV,,  . .  VB) 

and  the  functional  dependence  is  invertible. 


i  =  l,  ....,n 


(20) 


& 


ol 

$ 
I 


•S? 


Reconstruction  of  States  from  Measurements 


Similar  to  equation  (1)  the  system  equation  is 

V  =  f(V,U)  Vc 9?  ,  Ue&“ 

y  =  h(V)  y  (2i; 

Where  y  is  the  measurement  scalar  function,  differentiating  it  successively  with 
respect  to  time  would  resuit  in: 

.  _  -*h  ^h  rN 


Substitution  for  V  results  in 


y  =  - 


dv  dt 


f  =  Vh  •  f 


Similarly, 


y  =  -f  ) 


=  VLf(h)  •  f 


y(n)  =  V  L“_l(h)  •  f 

Since  y  is  a  measurement,  then  all  derivatives  can  be  considered  available.  Therefore 
the  problem  of  observability  of  the  system  would  be  reduced  to  the  conditions  under 
which  all  state  variables  are  determined  by  the  virtue  of  set  of  equations  (22). 


v 

Vh 

J 

VL,(h) 

yin) 

VLr'(h) 

•  f(V,U) 


In  order  that  vector  function  f(V,U)  be  determined,  again  the  inverse  of  observability 
matrix  derived  earlier  for  transformation  to  observable  form  must  exit,  namely 

Vh 

VL,(h) 


det  0  =  det 


vLrHh) 


Hence 


f(V,U)  =  [O]-1  •  Y  (25) 

This  condition  is  similar  to  the  first  condition  for  the  existence  of  transformation  from 
nonlinear  form  to  canonical  form.  Based  on  (25),  the  state  variables  can  be 
determined  from  the  measurement  vector  Y  if 


this  condition,  though  not  related  to  the  existence  of  any  transformation  of  system  of 
(21)  to  observable  form,  guarantees  the  existence  of  at  least  one  set  of  state  variables 
as  functions  of  measurements,  and  input(s). 

V  =  V  (  Y  ,  U) 


Observability  Analysis  of  MPD  Engine 

From  conservation  laws  for  species  in  plasma,  a  one-dimensional  global  model  can 
be  derived.  In  this  paper  the  application  of  observability  analysis  to  the  global  model 
is  presented.  However,  similar  steps  can  be  taken  to  achieve  more  laborious  task  of 
determination  of  observability  criteria  for  species  state  equations.  Those  global 
equations  are: 


continuity: 


momentum: 


energy: 


1 

11 

d\i 

+ 

<}P 

rJx 

= 

fi 

(27-a) 

<  >U 

+ 

RT  do 

oB 

‘ 

—  fj 

(27-b) 

dt 

u - 

<:bc 

P 

rbc 

P 

(E-uB) 

dT 

rJT 

+ 

RT 

c>u 

K 

O2  T 

cr 

.l 

~  f3 (27-c) 

dt 

U  /be 

Cv 

dx 

+  - 

pCy 

dx2 

E(E-uB) 

PCy 

It  is  assumed  that  shear  stress  and  energy  dissipation  due  to  shear  and  collision  along 
with  radiation  pressure  and  radiation  energy  are  negligible.  Also,  in  compliance  with 
species  set  of  equations  and  for  simplicity  of  having  one  control  variable,  the  applied 
magnetic  field  is  assumed  to  be  small  and  hence  with  negligible  effects,  i.e. 

B  ^  0 

Hence  the  control  variables  are 

U,  =  B(E-uB)  ^  0  ,  U2  =  E(E-uBj  =  E2  (28) 

Two  cases  for  available  measurements  are  considered.  In  the  first  case,  density  is 
assumed  the  only  measurement,  and  in  the  second  case  temperature  is  chosen. 


Case  I:  Density  Measurement 


0  0| 


P 

u 

T 


p(x,t) 


For  this  system  the  observability  matrix  can  be  determined  from  Vh,  VLf,  VL2. 


(29) 


Vh  =  VC  V  =  C  «  1,0,0] 

T  r-L  t  <9u  ,  ilP 

Lf  =  Vh  •  f  -  -  p—  +  u— 

OX  (JX 

I  Pi  pu 2 

VLf  =  -  u,  +  —  u  ,  pt  +  -  ,  0 


where  index  "i"  is  defined  as  an  abbreviation  for  "i"  times  differentiation  with  respect 
to  x. 


Lf2  =  VLf  •  f  =  Uj  +  —  u  (upj  +  put  | 

0 1  '  ' 


Pu2  RT 

+  Pi  +  uuj  H - Pi  +  RTt 


p2  Rpi  T, 

VL?  =  K,(  )  ,  K,(-)  1  +  R-1 

Pi  P  T1  jj 


Hence, 


p2  pu2 

O  =  —  Ui  +  — u  —  p,  4-  -  0 

Pi  u, 

pa,  Rp!  RT, 

K.(0  K,(-)  P!  +  —  —  + 

u,  p  T, 

The  first  condition  for  the  observability  of  the  system  is 

det  0^0 

This  condition  is  satisfied  if  both  of  the  following  conditions  are  met: 

I  i  I  du 

PiUj  +  pu2  =  —  p—  #0 
'  >  ch c  &x 


I\  '  w» 

P iTj  +  pT2  =  4r  P~r  *0  (32) 

The  second  condition  given  by  (26)  [and  in  a  more  restricted  form  by  ( 19) j  can  be 
determined  for  the  system.  However  the  implicit  function  theorem  which  results  in 
condition  (26)  for  finite  dimension  system,  will  be  modified  for  infinite  dimensional 
system.  In  equation  (25),  the  vector  function  in  the  left  hand  side  (i.e.  f(V,  U)),  for 
finite  dimensional  system  would  be  a  nonlinear  algebraic  vector  function,  whereas  in 
the  case  of  infinite  dimensional  system  it  will  be  a  nonlinear  ordinary  differential 
equation.  Therefore,  the  implicit  function  theorem  will  be  considered  in  the  space  of 
1-Jets  [34].  The  set  of  singular  points  of  the  equation: 
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F(x,  V(x),  =  0 


dV 


with  P  ss  - ,  in  three  dimensional  space  (x,  V,  P)  of  jets  are  those  for  which  the 

dx 

tangent  surface  to  F{x,  V,  P)  =  0  is  normal  to  (X,  V)  -  plane,  i.e.  —  =  0.  Applying 


this  criteria  to  the  sets  of  equation  (25)  would  result  in  the  fact  that  the  second 
condition  will  always  be  satisfied.  This  condition,  states  that  for  every  input  to  the 
system  there  exists  an  equilibrium  or  null  solution  for  the  system,  which  is  predictable 
for  physical  systems. 


Consider  the  Venn  diagrams  in  Figure  (2)  of  set  "I",  satisfying  the  condition 


and  "II"  satisfying  the  condition 


,  i 

-.be 


ij 

On 

p— — 

¥*  0 

<  hi 

0 

1 

£T 

#  0 

'hi 

P  dx  j 

:  system  equations  (27 

On 

_  _0_ 

c)T 

Ox 

=  0 


with  the  system  is  found  (i.e.  solved),  the  resulting  manifold,  i.e.  S  n  I  H  II  will 
represent  unobservable  manifold  of  S. 


From  equations  (31)  and  (32),  S  D  I  D  II  is  determined  by 


C'jl1)  n  r  dx 

Ui  =  -  .  <t>  X,t)  =1  — 

P  O  P 


u(x,t)  =  Cu<£(x,t)  +  u0(t) 

cT(t) 


Ti  = 


T(x,t)  =  CT(t)  0(x,t)  +  T0(t) 


(33) 


where  u0(t)  and  T0(t)  are  the  given  boundary  conditions  at  x  =  0  for  velocity  and 
temperature.  The  time  functions  Cu(t)  and  CT(t)  are  constants  of  integration  with 
respect  to  x.  Solution  of  (33),  by  considering  the  MPD  system  equations,  results  in  the 
following  observability  criteria: 


C  9  ^ 


C^(f>  +  +  T0 

+ 

Cx  R  Cu  K 

(CJ  +  u0)  +  “  ;  (CT<^  4-  T0)  -  ^  PlCT 

P  Cv  P  P°  Cy 

<r/p2Cy 


(34) 


However,  U2  is  the  control  variable  or  input  to  the  system.  Therefore,  to  satisfy 


observability  criteria  for  MPD  Engine,  using  density  as  the  measurement,  the  control 
input  should  satisfy  (34). 

Case  II:  Temperature  Measurement 

If  temperature  is  the  available  sensory  information,  then  the  measurement  function 
of  equation  (29)  will  be  modified  to 

P 

y  =  [0  0  Ij  u 
T 

the  same  steps  to  construct  the  observability  matrix  are  followed. 

0  0 

RT  u2 

(<tU2  -  KT2)  -  T,  +  -^7  - 

L  v  ul 

K,  Ku 

Rank  condition  requires  that 

K„  RT  u2 

-rr-  (®u.  -  kt,)  -  K,  t,  +  *r-  ~  *  0  <37) 

P  Cy  tyy  U, 

where  and  Ku  are  notations  used  for  very  complex  and  lengthy  terms. 

Due  to  the  complexity  of  the  condition  with  this  measurement  the  identification 
scheme  with  density  measurement  would  result  in  faster  and  less  complex  process 
which  makes  density  measurements  more  attractive. 

Conclusion 

The  criteria  for  observability  of  finite  and  infinite  dimensional  system  in  terms  of 
transformation  of  nonlinear  systems  to  observable  canonical  form  was  studied. 
Similarly  a  less  restrictive  approach  applicable  to  infinite  dimensional  systems,  which 
results  in  general  conditions  for  observability  of  such  systems  was  presented. 
Application  of  this  criteria  to  a  global  model  for  MPD  system  was  developed  and  the 
region  in  time-spatial  domain  for  control  action  was  derived.  Control  action  should 
avoid  this  region,  so  that  the  system  behavior  can  be  predicted  from  available 
measurements.  Selection  of  measurements  based  on  less  complex  criterion  to  satisfy 
was  presented  for  the  case  of  density  measurements  versus  temperature  measurements 
in  an  MPD  engine. 
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Figure  1.  Schematic  of  model  for  energy  balance  in  the  vicinity 

of  a  cathode  of  a  MPD  engine  in  the  absence  of  an 
applied  magnetic  field.  (1)  solid  conductor;  (2) 
conductor  in  "mushy"  state;  (3)  Debye  region;  (4) 
collisionless  region  outside  sheath;  (5)  Region  within  a 
recombination  length;  (6)  thermal  boundary  layer;  (7) 
free  stream. 
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Schematic  Representation  of  Observability  Region 
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NOMENCLATURE 

Description 
magnetic  field 

particle  velocity  in  laboratory  reference  frame 

peculiar  particle  velocity 

applied  electric  field 

body  force  on  particle 

current  flux 

heat  flux 

mas3  velocity 

diffusion  velocity  in  fluid  reference  frame 
specific  heat 
electron  charge 
energy 

Planck  constant 
Boltzmann  constant 
mass 

number  density 
pressure 

phase  change  energy 
radiation  energy 
temperature 
time 

electric  potential 
potential  energy 


X 


position 


rr 


density- 

electrical  conductivity 


mobility 


M 

r 

C. 

4 


mean  collision  frequency  of  type  s  particles  with  type  k  particles 

shear  stress 

particle  charge 

Debye  length 

mean  free  path 


Definitions 


b  - 


B 

B 


C,  =  mean  particle  speed  = 


8kT, 


|I/2 


-m. 


S'  =  S  +  Tf  X  B 


?=  S' 


nee 

f,  =  velocity  distribution  function 
P,  =  n,m,  =  species  density 
P  -  £  P, 


3 


m,mk 

ra.1.  =  reduced  mass  =  - 

m,  +  mk 


n,  =  generation  rate  of  species  s 


=  recombination  length  = 


D. 

D, 

De 


=  ambipolar  diffusion  coefficient  = 


T 

=  ion  diffusion  coefficient  =  k — 


Me  D,  +  Mj  Dg 
Me  +  M, 


=  electron  diffusion  coefficient  =  k — —  /< 
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IS, 


/ze  =  electron  mobility  = 


me  'eH 


•  Ma  ;/en  ^’i 

/Uj  =  ion  mobility  =  — — —  + 


P  ueH  P  P'« 


Min  - 


Mi.  = 


^ia  "ia 
e 


mi  ^eH 
us  =  u +0, 


t  —  mean  free  path  = 


DQ  «  ill 


Dt 


dt 


+  u  ■  V(  ) 


APPENDIX  I 


MODEL  OF  PLASMA  IN  THE  VICINITY  OF  AN  ELECTRODE 


A  partially  ionized  two-temperature  non-equilibrium  plasma  Sow  is  considered, 
which  obeys  Maxwell-Boltzmann  statistics.  The  induced  magnetic  field  is  included  for 
consideration. 

The  plasma  is  assumed  to  be  composed  of  three  monoatomic  ideal  gasses,  electron, 
ion  and  neutral.  The  pressure  of  each  gas  species  is  given  by  P,  =  n,kTs.  The  total 
gas  pressure  is  then  given  by 

p  -  £  Ps  (Li) 


The  energy  of  each  species  is  the  sum  of  the  translational  and  the  internal  energies. 
The  translational  energy  of  each  species  is  given  by 


3  P»kT, 
2  m. 


(1.2) 


A  Cartesian  orthogonal  coordinate  system  is  used  with  coordinates  that  are  parallel  to 
the  magnetic  field,  perpendicular  to  the  magnetic  field  and  perpendicular  to  the 
magnetic  and  the  electric  fields. 


The  electron  current  and  heat  flux  are  given  by  the  following  equations.  The 
electron  transport  properties  are  significantly  affected  by  the  magnetic  field  and 
therefore  these  effects  are  included. 


Je  =  0ii  en  +  fj.  +  crH  b  X  e  +  0 ii  V|,  Te  +  0±  Te 

+  Onbx  VTe 


(1.3) 


<le 


-  7  7'  -  v  1  v"  T>  -  x’*  T-  -  yR  K  *  VT- 

l  e 

—  Te  0,|  e ,  —  Te  Oj_  —  Te  Ojj  b  X  e 


(1.4) 


The  transport  coefficients  (<r,  o,  and  X')  are  presented,  for  example  in  Mitchner  and 
Krugerfl.l],  in  the  form  of  integral  functions  of  C.  The  transport  coefficients  are  also 
available  in  other  forms  from  sources  such  as  Bose  [1.2]. 

The  heavy  particle  transport  properties  are  not  affected  by  the  presence  of  a 
magnetic  field  unless  the  field  is  very  strong.  For  the  system  considered  the  magnetic 
field  is  assumed  to  be  sufficiently  weak  so  that  the  heavy  particle  transport  equations 
can  be  written  for  a  partially  ionized  plasma  without  a  magnetic  field.  Accordingly, 


the  following  relations  may  be  written  for  the  heavy  particle  transport  properties. 


du® 

c?u , 

i  .  -  — 

dxj 

dxn 

%  =  |  kTk  £  u,V,  -  V  VTb  -  ohkTh  £ 

£  a  s 


7h  =  Tj  -  en.Uj 


ff. 


2 

—  £  mkD,k  3k  -  —  D,t  V(lnTh) 

Q,P  k  /'* 


3,  >  V(i)  +  -  ^)V(lnP)  -  ^-(f-F,  -  V  nkFk) 

n  n  P  Ppm,  k 


where 

5n  y  =  Kronecker  delta  =  0  if  <y  ^  J  ,  1  if  a  =  J 

D3|t  =  concentration  diffusion  coefficient 
D,  =  thermal  diffusion  coefficient 


(1.5) 

(1.8) 

(1.7) 

(1.8) 

(I-®) 


Since  a  non-equilibrium  plasma  is  considered,  a  relationship  for  the  generation  of 
species  is  required.  Only  collisional  reactions  will  be  considered.  Specifically  the 
three-body  recombination  reaction  where  the  third  body  is  an  electron  may  be  written 
as  follows. 


e  +  A+  +  e 


e  +  A 


The  generation  rate  equation  is  then  given  by  the  following. 


ae  -  -T L  -  *(T.) 


(nen.)  -  nen,  \ 


(1. 10) 


(1-11) 


[•lia] 


illil] 


v  e  1/  u  I  i  '  7 

V  Qn  Jequil 

a(Tj  -  a,  1.09  X  1<T20  Te~9/2  (M3) 

=  Hinnov-Hirchberg  recombination  coefficient 


The  equilibrium  concentrations  are  given  by  the  Saha  equation,  namely. 


n«ni  1 

g, 

|  2Tm.kTe  | 

2 

2 

(  -ei 

nn  I 

=  2  — 

equii 

I  h2  | 

exp  ' 

[kT, 

(1.14) 


where  g,  is  equal  to  the  electron  energy  partition  function  and  ei(  the  ionizational 
energy.  Since  a  singly  ionized  gas  is  considered,  the  electron  and  ion  generation  rates 
are  equal  (ne  =  hj). 
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APPENDIX  II 

ENERGY  BALANCE  IN  THE  VICINTTY  OF  A  CATHODE 

The  energy  equations  are  written  specifically  for  each  region  shown  in  Figure  1.  In 
all  regions  the  plasma  is  assumed  to  be  non-neutral  (ne=/=nj)  and  two  temperature 
governed  (T,  ^  Tj). 

REGION  1  (solid) 


The  energy  equation  in  the  solid  material  is  given  by  the  well-known  Stephan  or 
heat  equation,  namely. 


E,t  -  (Eia  -  Eout)  +  Egen 


(II.  1  a) 


or 
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^solidC 


PhIM 


dt 


-  ^Ud^T)  + 


J2 

<7,oUd 


(Il.lb) 


where  Est  is  the  energy  stored  in  the  material,  (E,a  —  Eout),  net  energy  transfer 
through  the  material  by  conduction,  and  Es,n,  heat  generated  within  the  material  due 
to  Joule  heating. 

The  current  flux  is  given  by 

7  =  <7So, id^  (11*2) 


The  energy  input  to  the  solid  region  is  determined  from  an  energy  balance  at  the 
surface  y  =  0.  The  energy  balance  yields  namely, 


[\oiid 


<9T, 


solid 


dy 


=  [X 


<rr 


mushy" 


mushy 


<)y 


-o  +  H, 


i2L 

1  dt 


where  Hm  is  the  latent  heat  of  fusion  of  the  material 


(n.3) 


REGION  2  (Mushy) 


Region  2  is  the  "mushy"  region  where  solid,  liquid  and  gaseous  states  may  be 
present  simultaneously.  The  energy  for  the  "mushy"  region  is  similar  to  the  energy 
equation  for  the  solid  region  except  that  an  additional  terrn  is  needed  to  account  for 
the  energy  associated  with  the  material  phase  change,  Q  (T).  If  the  material  is 
assumed  to  "pure"  then  there  will  be  no  conduction  across  the  mushy  region.  The 
region  will  be  at  a  uniform  temperature  because  it  is  undergoing  a  phase  change.  If 
the  material  is  not  considered  "pure"  then  there  can  be  conduction  across  the  region 
because  there  may  be  a  small  temperature  gradient.  The  energy  transfer  in  the  mushy 
region  is  then  given  by  the  following  model  equation  (Reference  II.  1.),  which  is  in  the 
nature  of  a  modified  Stephan  equation. 

-  V(Xna!hyVT)  +  -2—  +  Q(T)  (II.4) 


where 


~  ^mushy^  (H-5) 

and  <7mushy  =  ^mushyl^)  =  electrical  conductivity  of  the  mushy  region  which  is 
temperature  dependent 
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The  two  surfaces  of  the  mushy  region  at  y  =  0  and  y  =  yt,  are  not  fixed.  The 
surface  at  y  =  0  is  allowed  to  recess  into  the  solid  material  as  melting  occurs.  Melting 
will  occur  in  the  region  as  long  as  the  energy  input  to  the  region  (at  the  surface  y  =  yt 
and  from  Joule  heating)  exceeds  the  energy  used  in  the  phase  transition  plus  the 
energy  removed  to  the  solid  region  by  conduction.  If  the  energy  inputs  and  outputs 
are  equal  the  boundary  will  not  recede.  The  surface  at  y  *  yj  will  recede  due  to  the 
evaporation  of  material  at  the  surface. 

The  energy  input  to  the  mushy  region  from  the  plasma  is  determined  from  an 
energy  balance  evaluated  at  the  surface  y  =  yt,  as  shown  in  Figure  1.  The  cathode  is 
assumed  to  have  a  catalytic  surface  where  incident  ions  and  electrons  recombine  and 
are  re-emitted  as  neutrals.  Thermionic  emission  may  also  occur  if  the  cathode 
temperature  is  sufficiently  high.  It  should  be  noted  that  the  thermionically  emitted 
electrons  provide  an  additional  localized  current  which  in  turn  produces  an  additional 
localized  Joule  heating.  The  localized  heating  may  lead  to  localized  evaporation  or  to 
the  eruption  of  material.  The  overall  energy  balance  equation  is  then  as  follows. 


Q  "h  dcond  dp  "b  dr*d  "b  Qj 


(n.6) 


where  qp  is  the  net  energy  transfer  to  the  surface  from  the  plasma  due  to  the  kinetic 
energy  and  potential  energy  of  the  particles,  qr4d,  the  net  radiation  to  the  surface  from 
the  plasma,  Qsut,  the  energy  generation  at  the  surface  due  to  recombination  of  ions 
and  electrons  Q  ,  the  energy  associated  with  phase  change  and,  4:ond*  energy 
conducted  into  the  mushy  region. 


REGION  3  (Sheath) 

The  region  immediately  adjacent  to  the  electrode  and  within  a  distance  of  the 
order  of  Debye  length  (<p)  of  the  surface  is  the  sheath  region.  In  this  region,  charge 
separation  occurs  and  a  net  negative  charge  exists  because  of  an  excess  number  of 
electrons.  The  region  is  considered  collisionless  in  the  sense  that  only  electron-neutral 
and  ion-neutral  collisions  are  present.  These  collisions  are  included  because  of  the 
large  number  of  neutrals  in  the  region.  Since  neutrals  being  emitted  from  the  surface 
do  not  experience  a  force  from  the  fields,  they  tend  to  remain  near  the  surface.  They 
are  removed  through  diffusion  driven  by  the  concentration  gradient.  All  other 
collisions  are  assumed  to  be  negligible. 

The  energy  transfer  to  the  surface  from  the  plasma  is  from  the  impact  of  particles 
on  the  surface.  Since  a  cathode  is  being  considered  the  particles  of  interest  are  the 
ions.  The  energy  of  the  particles  is  in  two  forms,  the  kinetic  energy  due  to  their 
motion  and  the  potential  energy  associated  with  moving  charged  particles  through  an 
electric  potential.  The  electric  potential  will  tend  to  move  the  electrons  away  from  the 
electrode  while  accelerating  the  ions  towards  the  electrode.  The  neutral  particles 
possess. 

Since  the  interest  is  in  particles  which  strike  the  surface,  the  velocity  component 
normal  to  the  surface  is  the  one  of  interest.  The  describing  energy  equations  are  given 
by  the  following. 
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The  species  heat-flux  vector.  qfs,  is 


%  -  /_“  J  n,  m,  C2  Cx  f,  d3 


(H.7) 


and  the  particle  potential  energy,  Ws,  is 


W. 


-  -v*dy 


in.*) 


where  the  electric  potential  V,  is  given  by  the  following  Poisson's  equation. 


VJV,  -lx  f,f,d3c 

The  species  current-flux  vector,  T,,  given  by 


(U.9) 


—  J_x  Ds  fs  f,  ^ 


(II- io) 


The  distribution  fs  is  determined  from  the  Boltzmann  equation,  namely 


REGION  4 


5  I 


Du 


<*f. 


“(Qsfs)  +  (u/  +  Kf*)  + 

~  c’"’  TcT  “  v""1 


m,  Dt 

« T,  -?u. 


*  :>c 

= 


(ii-ii) 


where  is  the  rate  of  increase  of  the  property  of  interest  (mass,  momentum,  energy 
or  charge)  due  to  collisions  between  particles  of  type  s  with  particles  of  type  k 


Region  4  is  very  similar  to  region  3,  as  they  are  both  considered  collisionless.  The 
describing  equations  are  therefore  the  same  for  the  two  regions.  Since  this  region  is 
outside  of  the  sheath,  it  is  expected  to  contain  a  greater  number  of  ions  than  region  3. 
The  electron-neutral  and  ion-neutral  collisions  also  become  negligible  because  of  a 
large  decrease  in  the  number  of  neutrals  in  the  region.  The  collision  parameter  4^  in 
equation  II.  11  is  therefore  equal  to  zero. 


REGIONS  5  and  6 
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Regions  5  and  8  are  assumed  to  be  collision-dominated  and  can  therefore  be 
described  using  the  hydrodynamic  approximation.  Particle  recombination  is  assumed 
to  absent  in  region  5  because  it  is  within  a  distance  of  the  order  of  the  recombination 
length  2r  of  the  electrode  surface.  The  recombination  is  included  in  region  6.  The 
resulting  energy  balance  equations  for  Region  5  are  as  follows. 


~(y  n,kT,)  +  (f  nakTe)V-u.  -  -Vq.  +  T,:VZt  +  7.-EP. 


2  m.  3 

-  —  7  k(T.  -  Th)  -  R. 

mh  2 


(11.12) 


^-(fa„kTb)  +  inkkTbV-irb  -  +  rb:Virb  +7-f, 


2m.  3 

-  — 2-  uhenhf  k(Th  -  T.)  -  Rb 
mb  2 


(n.i3) 


where 


nVu  -  r „  ;  nh  -  n,  +  nu  ;  mn  m,  =  mh 

'  TC  / 


The  transport  equations  are  given  in  appendix  I. 


The  energy  equations  for  region  6  are  given  by  the  following.  The  recombination 
energy  is  included  in  the  electron  energy  equation. 


kT.  +  «,))  +  Q.(|  kT,  +  €,)V-U,  -  -V-ra,  -  ^7.) 


-»  _  2m.  3 

+  V«i  +  - 1  7,„n,  fklT.-TJ 

mh  2 


(11.14) 


nhkTh)  +  y  nhkTbV  uh  =  -V  qh  +  \.Vnh  +  7.-E', 


-  — =■  ~h«nh  7  MTh  -  T.)  -  Rh 

mh  2 


(11.15) 


v'.aV'j.  -  j  j  .•  .-V  j.  j\r  .■  V'-V.  v 


\  *v  s 


REGION  7  (Free  Stream) 


In  the  free  stream  region,  the  various  gradients  are  assumed  to  very  small 
compared  to  the  other  terms  and  are  neglected.  The  rate  of  change  of  energy  within 
the  system  is  equal  to  the  generated  energy  from  Joule  heating  minus  the  energy  lost 
from  radiation,  recombination  and  collisions.  The  energy  equations  can  therefore  be 
simplified  as  follows. 

,o  o  2  m 

■£■(«.(  7  kT.  +  6,))  - 1  uAn.  7  k(T,  -  Th)  -  R.  (11.18) 

ca  l  mb  l 

77 (j  nhkTb)  -  —  Phenh  ±  k(Th  -  T.)  -  Rh  (11.17) 

7,  =  rrn  e„  +  e±  +-THbX(  (11.18) 
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ABSTRACT 

A  method  for  analyzing  controllability  of  non¬ 
linear  distributed  parameter  systems  is  presented.  This 
analysis  is  based  on  the  Lie  algebra  augmented  for  dis¬ 
tributed  systems.  The  augmentation  is  performed  by 
transformation  of  the  state  equations  into  a  finite 
dimensional  state  space.  The  resulting  controllability 
criteria  is  applied  to  a  nonlinear,  distributed  electric 
propulsion  system,  namely,  a  magneto-plasma  dynamic 
(MPD)  engine. 

INTRODUCTION 

Distributed  parameter  systems  are  often  approxi¬ 
mated  or  modelled  by  a  lumped  parameter  approach  or 
by  discretization  of  the  system  in  time  and/or  space 
[  1 .2) .  This  approach,  although  it  has  been  useful  for 
design  of  some  control  systems,  in  general  cannot  be 
applied  to  an  accurate  dynamic  analysis  of  nonlinear 
distributed  systems.  This  is  especially  true  when  con¬ 
siderations  of  general  stability,  controllability  or  obser¬ 
vability  of  the  system  is  concerned.  Therefore,  more 
accurate  methods  should  be  developed  to  investigate 
the  general  behavior  of  the  nonlinear  distributed  sys¬ 
tems.  Regarding  the  questions  of  stability,  many 
authors  have  studied  the  nonlinear  distributed  parame¬ 
ter  systems.  One  of  the  first  attempts  using  Lyapunov’s 
direct  method  for  infinite  dimensional  systems  was 
made  by  Massera  3  .  Based  on  an  abstract  theory  of 
Lyapunov  stability  scheme  of  infinite  dimensional  sys¬ 
tem.  references  4,5  have  defined  a  Lyapunov  function 
on  Hilbert  space  and  Banach  space  cases. 

For  finite  d  imensional  systems  the  concept  of  con¬ 
trollability.  was  introduced  in  early  1960  s  by  Kalinan 
and  others  for  linear  systems,  and  in  the  early  1970's  by 
work  of  Herman  6  and  Haynes  and  Hermes  7  for  non¬ 
linear  systems.  The  nonlinear  analogy  of  linear  control¬ 
lability  criteria  for  lumped  parameter  systems  is  carried 
out  independently  by  Sussman-Jurdjevic  81  and  Kroner 
9  by  applying  the  Lie  Algebra. 


This  study  is  supported  by  the  f.S.  Air  lores  Office  ..(  Srin.liGr 
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Controllability  of  a  linear  distributed  parameter 
system  has  been  investigated  by  researchers  10,1 1,12] 
for  the  cases  that  the  control  energy  is  either  distri¬ 
buted  in  the  system  or  is  present  at  the  boundary.  As  a 
recent  attempt  to  derivation  of  an  exact  controllability 
criteria,  Fattorni  and  Russel  1 1 3 1  have  used  reduction  of 
a  special  parabolic  system  to  a  moment  problem  and 
related  the  controllability  to  the  absolute  convergence 
of  the  moment  problem.  • 

In  this  paper  a  transformation  is  presented  to 
augment  the  state  space  of  nonlinear  distributed 
parameter  system.  This  transformation  reduces  the 
system  into  an  equivalent  set  of  finite  dimensional  non¬ 
linear  systems.  The  method  proposed  by  [14,151  which 
transforms  the  finite  dimensional  nonlinear  system  into 
a  canonical  form  is  employed.  Based  on  this  transfor¬ 
mation,  a  controllability  criteria  for  nonlinear  distri¬ 
buted  parameter  systems  is  established. 

SYSTEM  REPRESENTATION 

Distributed  parameter  systems  can  be  represented 
in  many  different  forms.  However,  in  terms  of  energetic 
systems  where  conservation  laws  and  Maxwell  equations 
may  be  applied,  a  distributed  parameter  system  can  be 
represented  by  the  following  form  of  state  equations. 

~  =  A(V)  +  B( V)  U  (1) 

where  V(x,l)  is  a  state  vector  in  a  Banach  space,  which 
is  a  function  of  time  (t)  and  space  dimension  (x),  and  it 
is  composed  of: 

V(x.t)  =  Col  [ V i ( x , t ),  v2(x.t) . Vfl(x,t)]  (2) 

The  operator  A  is  a  nonlinear  operator  defined  on  the 
Banach  space.  The  specific  case  considered  here  is 
when  A  forms  a  general  nonlinear  function  F*  of  the 
following  form: 

.  'V  ,  <mV 

A ( V )  =  F‘  V,  — .  —  (3) 

,  >x  .  'Xm 
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The  input  space  consists  of 


L'  =  col 


<-  i(x.t)  .  Uj(x.t) .  Ur(x,t) 


(«) 


The  augmented  state  and  input  space  are  defined  as 

f  r  I 


Therefore,  the  operator  B(V)  will  be 


V 


B(V) 


B,(V)  B.(V)| 


Br(V) 


(5) 


i  id) 


with  B,(V),  1  <  i  <  r,  being  a  nonlinear  operator  of  the 
form 


B,( V)  =  g/  V. 


y_ 


■p\' 

*xp 


(6) 


V*"11 


'mr 


with  p  <  in.  Equivalently,  the  general  operator  B(V) 
can  be  represented  as 


B(v)  =  G* 


■  •py 

-*xp 


(7) 


where 


G 


Equations  (1),  (3)  and  (8)  represent  a  wide  class  of  non¬ 
linear  distributed  parameter  systems.  This  system 
classification  and  representation  is  used  in  the  following 
controllability  derivations. 

SYSTEM  AUGMENTATION 

In  order  to  derive  controllability  criteria  for  the 
system  of  equation  (1),  the  system  is  augmented  such 
that  the  resulting  equivalent  system  can  be  formulated 
in  terms  of  a  finite  dimensional  nonlinear  system.  This 
augmentation  will  reduce  the  infinite  dimensional 
domains  of  the  A  and  B  operators  into  a  finite  domain 
at  each  point  in  the  space  dimension  x. 

Let 

—  -V»>,  =  V(2»  .  •  -  •  (9) 

he  be2 

Therefore, 

~  =  F*  |v,  V").  V'21 . V,m) )  (10) 

+  G*  (v,  V»>.  V12’ . V'PI  jy 


Differentiation  with  respect  to  x  of  the  above  state 
equation  results  in 


“x 


=  —  p*  y  vll(  V1"’1 
■  *x  1 . 


G*  \.  Vu) . v,pl  t; 


HD 


Assuming  V,  F*.  and  G‘  are  analytical  functions,  i  hen 
equation  (11)  can  be  written  as: 

,  i\’"l  -F> 

-  =  --  (')  +  — -  U  +  G*  — 

■  't  ‘X  "X  'X 


Differentiating  the  state  equation  m  times  resuits  in 
.  <V<2>  .<*F*  ,  ,2G  ’  ,.  2  ‘G  *  r  G-  '21. 

—  -  v"<,  +  “ 1  f  — -y - " 

,  (yto'l 
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fill 
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Or  in  the  augmented  space  the  system  of  equation  (11) 
can  be  written  as 


T  =  F(M  +  c: 


115) 


where 


F(  ‘ )  = 


F*(M 

r  f*(  ') 


(16) 
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The  augmented  form  for  G*.  G.  would  he 
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Equation  (151  is  an  augmented  form  of  equation  (l|  and 
(101  which  can  be  represented  in  an  augmented  lini'e 
dimensional  space  at  each  x  with  F:  x  —  \  and 

DERIVATION  OF  CONTROLL.VB1LITY  CRITERIA 

In  the  previous  section,  it  was  shown  that  a  non¬ 
linear  distributed  system  can  be  represented  in  an 
equivalent  nonlinear  finite  dimensional  system  in  the 
following  form. 

:(t)  =  F I,  I  —  til  Ol  (tl  1 18) 


where  '  £  R4  .  I  €  R  q  =  nm 

In  general,  controllability  concept  is  related  to 
relationship  between  inputs  and  slates  16  .  Namely, 
each  state  variable  ■  £  M  can  be  manipulated  by  an 
input  I  £  11.  In  order  to  establish  this  relationship,  the 
system  equations  should  be  transferred  into  a  canonical 
form,  as  shown  in  the  following. 

=  f!  •*)  -t-  gl  (19) 


and  g  = 


The  transformation  that  generates  ■*  •  M*  from  ■  •  M 
is  a  one-to-one  mapping  T|  •■):  M  — >M ‘ .  Hence 


'*  =  T  | ' ) 


N1  C  Rq 


and  inverse  mapping  is  denoted  by  T(  v*):.Vl*— *M 

■  =?('*)  ‘  •  M*  C  R4  (21) 

Therefore. 

“t  =  ~r~  [f  +  G1J22) 

Comparison  of  equations  (19).  (20)  and  (22)  results  in 

T.l  •  | 

T  1  si  '  I 

—  F=flg-|=  (23) 


Since  T,(  ')  is  a  scalar  function,  then 

T  1  -T 

- -  p  =  \'  -  F  =  'CT.F  =  vdT,.F>  (25) 


Therefore 


L.,  =  < dT  .  F>  for  i  =  l . q-1 


T,  =  <dT,.  F > 


Tj  =  <  d<dT, .  F  >.  F  > 


Furthermore,  combination  of  equations  (19)  and  (22) 
vields 


- C(  )  =  g  = 


or  equivalently 

<dT,.  G(  ')>  =  0. 


i  =  1 . q-1 


Combination  of  (26)  and  (27)  results  in 

<dTl_1,  C(  ')>  =  <d<dT„  F>.C> 

By  direct  expansion  it  can  be  proved  that 

<d<dT,.  F>,  G  >  =  —  <dT,.  F  -  C 
4-  <d<dT,.  G>,  F> 


'T  '  'F 

The  vector  -  F  —  -  G  is  a  Lie  bracket  and  can 

. t\  . 1  \ 

equivalently  be  represented  by  the  Lie  bracket  F.G  or 
(ad  sup  1  F.G).  Moreover,  it  was  shown  by  (27)  that 
the  inner  product  <dT,.  G>  is  zero.  Therefore. 


<dTj,  G(  x)>  =  <d<dT,.  F|  ■)>, 
G(  ■)>  =  -  <dT,.  ad'F.GO  =  0 


Similarly. 


<dT..,,  G(  'l>  =  <dT  .  ad  F.G 


Since  all  of  the  inner  product  terms  have  the  common 
'T, 

vector  dT,  =  — — ,  then  the  above  results  can  be  sum- 

mariied  in  a  matrix  notation  as  follows. 

4T  i 

G(  Ml(ad1F,G)|(ad'’l',(.;j:  •••  |(adq  *F.G) 

=  0,  0 _ I  — 1  )q  1  (28) 

Let  define 

C  =  G  I  ( ad*F.G )  j  •  •  •  ladq_IF.G)  (29) 


The  Momentum  equation  can  be  written  in  the 
x-direction  as 

,  .  #ii  #P 

,, -H.  +  ,,u  =  -  -  +  JB  (31) 

,  >t  #x  *x 


It  is  assumed  that  — ^  =  0,  and  the 

,  #v  'z 


electromagnetic  force  J  X  B  is  applied  in  x  direc¬ 
tion.  J  is  the  local  current  between  the  elec¬ 
trodes. 


J  -  J(x,t) 


Then  (281  can  be  written  as 


,#T 

#T 

0.  0.  . 

.  l-l)q‘‘ 

C"‘ 

f cV  “T  +  u 

dt 

•  't 

In  order  to  have  controllability  for  the  system 

*T, 

represented  by  (18)  or  (19),  -  has  to  exist.  This 

requires  that  the  controllability  matrix  C  has  to  be 
invertible.  Therefore,  controllability  analysis  of  the 
nonlinear  distributed  parameter  system  would  result  in 
checking  the  rank  of  the  C  matrix  shown  by  equation 
(29). 

APPLICATION  TO  MPD  ENGINE 

t  -  System  Formulation 

Recent  increase  in  space  missions  and  construction 
of  the  space  station  has  attracted  attention  to  new 
alternatives  to  chemical  propulsion  systems.  One  such 
system  is  an  electric  propulsion  engine  or  magneto¬ 
plasma  dynamic  (MPD)  engine  17 

Consider  a  simplified  model  for  a  Magneto  Plasma 
Dynamic  Thruster  which  has  the  applied  magnetic  field, 
B.  and  applied  electric  field.  E,  as  the  input  control 
variables.  The  general  governing  conservation  equa¬ 
tions  for  the  system  can  be  written  as 


Continuity  equation 
D,  ___ 
Dt  1 


~  ,T  -  ii  «  0 


f  or  a  one-dimensional  model  the  above  equilion 
can  be  reduced  io: 


-  -  —  i ,  u )  =0 

't  dt 


One-dimensional  non-steady  energy  equation  wil 
lead  to: 


=  -P  —  +  JE 

■  dt 


In  general,  pressure  P  is  a  nonlinear  func¬ 
tion  of  density,  /#,  and  temperature,  T,  i.e.  the 
state  variables.  In  case  of  perfect  gases,  this 
relation  is 

p  =  P(,..T)  =  ,'RT  (33) 

In  the  energy  equation  it  is  assumed  that 
diffusion  term  is  negligible  compared  to  the  heat 
generated  inside  the  flow  and  the  energy  resulted 
from  compressibility  effects  18j.  The  applied 
fields  are  related  to  the  current  density  by  the 
Ohm's  law; 

J  =  .ifE  —  uB)  (34) 

where  the  parameter  "r  is  electrical  conductivity 
and  it  is  a  property  of  the  plasma. 

Rewriting  equations  (30),  (31)  and  (32)  with 

respect  to  substitutions  of  equations  (33)  and  (34) 
resuits  in  a  set  of  state  equations  as: 

—  +  —  (/'u)  =  0  1 35a) 


'U  *u 

-  4-  U  - 

't  dc 


R_  '( , ;T )  _  -Dili  -  uBI 
r  *x  = 


'T  'T  RT  'u  -EIE  -  uB)  . 

-  +■  u  -  -i-  -  -  =  - ‘ - -  1  3.ic  I 

'*  ■  dt  .  cv 

The  set  of  equation  135)  can  be  formulated  in  a  distri¬ 
buted  space  form  as: 


where  ,  and  u  are  the  plasma  density  and  velo¬ 
city.  respectively. 


-  G*|Vji; 
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The  control  action  U  for  this  formulation  will  be 

B(E-uB)]  [u, 

U  = 

E(E  —  uB)  [UJ. 

The  distributed  state  variables  are 


V  =  u 
T 


11  -  Controllability  Analysts 

For  the  system  of  equations  (35)  the  controllabil¬ 
ity  technique  presented  in  this  paper  is  applied.  In 
relation  to  the  form  of  equation  (36).  it  can  be  shown: 
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Therefore,  from  (12),  (13),  (16)  and  (17)  it  can  be  shown 

that 

Based  on  equation  (29).  the  controllability  matrix  C  can 
be  constructed. 

C  =  I  g2  j  g3 1 |  {ad'F.gi )  j  (ad'F.gj)  j  (ad'F.gj)  j 
(ad’F,g4)  j (ad2F.g[)  j  (adJF.g2)  j  (adzF,g3)  J  (adJF.g<) 


The  rank  condition  for  this  controllability  matrix 
would  lead  to  the  criteria  for  the  system  controllability. 
As  shown  before,  the  terms  in  the  controllability  matrix 
(39)  can  be  calculated  as: 
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The  same  would  be  considered  for  glt  g3  and  g4  in  the 
calculation  of  ad1  and  ad2. 
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Following  the  detailed  calculation  to  establish  the 
matrix  C,  the  rank  condition  of  matrix  C  would  be 
satisfied  if  and  only  if 

rr8  An  ,  l 

-L.  —  (pu)  *  0  (40) 

p8  Cy  'Ot  'AX 

Based  on  the  fact  that  cr  1S  not  zero,  both  #  0 

'>X 

and  (pu)  #  0  provide  the  controllability  conditions 

for  the  operation  of  the  MPD  engine.  Therefore,  to 
have  a  controllable  engine  neither  density  nor  momen¬ 
tum  of  the  fiow  should  be  constant  along  the  engine 
axes. 


CONCLUSION 

Controllability  of  nonlinear  distributed  parameter 
systems  was  studied.  A  criteria  for  controllability 
analysis  was  developed.  This  technique  is  based  on 
applications  of  Lie  algebra  for  augmented  nonlinear  dis¬ 
tributed  systems.  The  augmentation  method,  which 
involves  two  transformations  of  the  distributed  equa¬ 
tions,  was  presented.  The  resulting  controllability  cri¬ 
teria  was  applied  to  a  magneto-plasma  dynamic  (MPD) 
engine.  This  electric  propulsion  engine  has  nonlinear 
characteristics  and  is  a  distributed  parameter  system. 
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Abstract 

An  electro-magnetor-plasma  process  involves  nonlinear,  coupled,  and  complex 
dynamics  and  is  subjected  to  a  number  of  instability  modes,  some  of  them  being 
associated  with  ionization  or  electro-thermal  instability.  Systems  involved  in  such 
processes  are  more  complicated  due  to  the  spatial  and  time  variations  in  the  system 
states  (distributed  parameter  system).  The  objective  of  this  paper  is  to  present  a  new 
approach  to  stability  of  such  distributed  parameter  systems.  This  technique  is  based 
on  spatial  extension  of  Lyapunov  stability  theorem.  This  Lyapunov  functional 
approach  is  applied  to  a  magneto-gas  dynamic  problem,  and  stability  results  are 
presented  and  discussed. 

Introduction 

Dynamics  of  systems  and  processes  which  involve  multi-energy  modes  and 
interactions  between  fluid,  thermal,  and  electromagnetic  fields  can  generally  be 
expressed  by  partial  differential  equations.  Electric  propulsion  systems  such  as 
Magneto-Plasma-Dynamic  (MPD)  engines  are  examples  of  such  systems.  Stability 
analysis  of  such  systems  is  very  complicated  by  the  fact  that  state  variables  are 
functions  of  both  space  and  time.  In  this  paper  a  new  approach  for  stability  analysis 
of  such  distributed  parameter  system  is  presented.  This  approach  is  based  on  the 
Lyapunov's  direct  stability  method.  Lyapunov  stability  theorem  has  become  an 
important  vehicle  in  derivation  of  stability  analysis  of  solutions  to  linear  and  nonlinear 
ordinary  differential  equations.  This  approach  attempts  to  make  statements  about 
stability  of  motions  of  a  dynamical  system  without  any  knowledge  of  the  solutions  to 
its  governing  equations. 

Although  the  development  of  Lyapunov’s  stability  theorem  for  ordinary  differential 
equations  has  been  widely  investigated,  its  application  to  solutions  of  partial 
differential  equations  (distributed  parameter  systems)  has  been  limited.  Most  of  the 
stability  results  for  distributed  parameter  systems  have  been  derived  by  use  of 
approximation  methods.  These  methods,  in  general,  use  reduction  of  the  partial 
differential  equations  to  a  system  of  ordinary  differential  equations  by  spatial 
discretization  or  by  assuming  a  harmonic  time  dependence.  In  the  harmonic  case  the 
Galerkin  method  based  on  a  truncation  of  the  modal  expansion  is  used.  There  have 

*  This  research  as  been  supported  under  AFOSR  Grant  86-0278,  with  Dr.  J.  Tishkoff  as  Technical 

Monitor  of  the  Project. 
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been  studies  in  the  applications  of  Lyapunov’s  stability  theorem  to  distributed 
parameter  systems  (DPS)  without  any  approximation.  First  attempt  to  apply 
Lyapunov’s  direct  method  for  DPS  was  made  by  Massera  [l)  and  Zubov  (2 j .  Massera 
extended  Lyapunov’s  method  to  denumerably  infinite  system  of  ordinary  differential 
equations  to  arrive  at  a  functional  in  terms  of  infinite  dimensional  state  vector.  This 
functional  was  used  in  place  of  the  Lyapunov  function.  There  have  been  studies  in 
formulation  of  Lynpunov  functions  for  special  problems  [3,4,5).  A  general  stability 
theorem  based  on  the  existence  of  a  Lyapunov  functional  is  established  by  Zubov  in  a 
functional  space.  Zubov  has  considered  the  general  type  of  system,  namely, 

^(t,X)  =  LZ(t.X) 

where  JJ(t,X)  is  an  n-dimensional  vector  valued  function  defined  over  some  region  17  of 
spatial  domain.  Matrix  valued  operator  L  is  a  differentia!  operation  defined  on  the  17. 

The  abstract  theory  of  Lypapunov  stability  of  infinite  dimensional  system  was 
studied  by  Buis,  Vogt,  Eisen  .6),  Pau  |7] ,  and  Banks  [8).  In  these  studies  Lyapunov 
function  is  defined  by  <  X,SX>  where  S  is  a  bounded  positive  self-adjoint  operator, 
for  Hilbert  state  space,  and  as  [  X.SXj  for  the  Banach  space  case,  where  [-,*]  is  a  semi- 
inner  product. 

Another  approach  to  the  stability  analysis  has  been  the  application  of  the 
frequency  domain  methods.  In  this  approach  a  generalized  circle  criterion  is  used  for 
infinite  dimensional  systems  [9, 10,11). 

I 

Lyapunov  Functional  Approach 

A  distributed  parameter  system  can  be  defined  by  the  following  equation 

'  ;i\  ,iDx 

|  -  J<  x.  .f.‘>  w 

)  where  Xis  the  state  vector,  a  is  spatial  coordinate,  Uis  the  state  input  and  t  denotes 

time.  If  a  Lyapunov  functional  — ►  di  is  defined  for  this  system  such  that  it 

;  satisfies  the  following 

i 

i.  V(0)  =  0,  and  V(X)  >  0  for  X=*=  0  and  Xeffi 

ii.  Vec^Sft0) 

iii.  V  =  <gradV,  _f  >  <0, 

then  the  following  two  theorems  can  be  stated  12:. 

Theorem  1:  If  in  a  neighborhood  of  the  state  origin  there  exists  such  a  Lyapunov 
functional  for  the  system  (I),  satisfying  (i)  through  (iii),  then  the  origin  is 
a  stable  state  for  the  system. 

Theorem  2:  If  in  a  neighborhood  of  the  origin  there  exists  a  Lyapunov  functional  V 
for  the  system  (1),  satisfying  (i)  through  (iii)  everywhere  except  at  the 
origin  itself,  then  the  origin  is  an  asymptotically  stable  state  for  the 
system. 
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Therefore,  based  on  these  theorems,  stability  of  distributed  parameter  systems  can 
be  guaranteed  provided  a  Lyapunov  functional  can  be  derived.  In  general,  for  systems 
defined  in  a  Banach  space,  the  functional  V  can  be  constructed  based  on  a  semi-inner 
product,  namely, 

V  =  [  X,  PXj 

and  for  systems  defined  in  a  Hilbert  space,  the  functional  can  be  formulated  by  inner 
product,  namely 

v  =  <  X,  PX>  . 

In  the  following  section,  derivation  of  such  functional  is  shown  for  a  simplified  model 
of  an  MPD  engine. 

Problem  Formulation: 

Consider  the  MPD  engine  shown  in  Figure  1.  The  plasma  dynamic  equations  for 
this  engine  consist  of  Maxwell’s  equations,  Ohm’s  law,  conservation  of  electric  charge, 
equation  of  state  (ideal  gas  law)  and  a  set  of  mass,  momentum  and  energy  equations 
[24 j .  These  equations  can  be  written  as  follows: 


Maxwell  equations: 


V  X  H  =J  + 


V  X  E  =  — 


V  H  =  0 


V  -E  =  -Pe 


Ohm's  law: 


J,  =  <j  j^E,  +  pe(U  X  H);  +  peU,  i  =  x,y,z 


Conservation  of  electric  charge: 


>p,  3  /»J 

-  +  v  — 1 

J-l  -*J 


Equation  of  state  (ideal  gas  law): 


P  =  RpT 


m 

9RR| 


Xi 

§8 


Conservation  of  Mass: 


^  4-  V  —  (pU:)  =  0 
j-i  '*j  J 


Conservation  of  Momentum: 


HP  3  Crj: 

-  —  +  v;  -21  +  FY  +  F. 


where, 


Fe.  =  pe  E,  +  /ie(J  X  H), 


F_.  =  gravity  force  per  volume 


Tn  = 


i  =  x,y,z  or  x1,x2,x3 


Energy  Equation: 


;,Pe m  N3,  ;iPem  Uj 

,:,t  ^  ;/xj 


'U.r,.  CQ. 

+  EjJj  + 

'  *Xj  J  i  hij 


where  em  =  total  energy  per  unit  mass  =  cvT. 

In  this  model  it  is  assumed  that  the  plasma  is  originally  at  rest  with  pressure  P0, 
temperature  T0,  and  density  pQ.  An  external  uniform  magnetic  field  H0  is  applied  to 
the  system  where, 

H0  =  i  Hx  4-  j  Hy  4-  k  0  . 

There  is  no  electric  field  applied  to  the  system.  Plasma  is  perturbed  by  a  small 
disturbance  and  as  a  result  the  state  of  the  system  is  a  combination  of  stationary 
(equilibrium)  part  and  perturbed  portion.  For  velocity  vector  the  basic  flow  is  zero. 
Therefore, 

U  =  i  u(x,t)  4-  jv(x,t)  4-  k  w(x,t) 

However,  an  assumption  is  made  that  the  variations  of  variables  are  only  functions  of 
one  spatial  dimension,  x,  and  time.  Therefore,  instantaneous  pressure,  temperature 
and  density  can  be  written  as 


P  =  P0  +  P'(x.t) 
T  =  T0  +  T'(x,t) 

p  =  p0  +  p'M 


>  , 
V, 


I 


t? 


V 

j 

i 


|C> 


j 


Electric  and  magnetic  fields  can  be  represented  as 

E  =  i  Ex(x,t)  +  j  Ey(x,t)  +  k  EJx.t) 
H  =  H0  +  h  (x.t) 


Hx  +  hx(x,t) 


+  J 


Hy  +  hy(x,t) 


4-  k  h,(x,t) 


current  density  J  and  net  electric  charge  pe  are 
J  =  J  (x,t),  pe  =  pe(x, t)  . 


The  one-dimensional  assumption  results  in:  =  0,  ( [  ^  —  0, 

(ry  (j  z 


Based  on  the  aforementioned  treatment  of  the  problem,  the  following  describing 
equations  can  be  derived. 


Maxwell’s  equations: 


<’>E 


Jx  +  e— JL 
x  <'t 

=  0 

(2) 

/'Ev 

4h, 

4-  c  — 

L 

(3) 

■  4 

'hi 

4E, 

i>  hv 

I  +  = 

1  <>t 

y 

/be 

(4) 

.'h„ 

A 

'V't  ~ 

0 

(5) 

4hv  ( 

'E, 

-  - 

(6) 

.  a 

'he 

'  ;h7 

i  'Ey 

z 

Pe  ~ 

y 

(7) 

'he 

Generalized  Ohm’s  law: 


Jx  -  <r  (Ex  -  wHy)  +  pe u 


(8) 


vj 


Jy  =  a  (Ey  +  Me  wHJ  +  pev 


Jt  =  <t  (Es  +  PeuH  -  Mev Hx)  +  pe w 


Conservation  of  electric  charge: 


'  Pe  '  'Jx 

— -  +  — 1  -  0 
<  't  <  *x 


Equation  of  state  for  perturbed  variables  is 


pi  o'  T' 

—  =  —  +  —  .  where  PQ  =  p0RT0 

Po  Po  T0 


Linearized  continuity  equation  becomes: 


tip1  .'>n 

~h~  +  Po~ =  0 

■  it  -be 


Linearized  equations  of  momentum  are: 


«"'P'  ,  4  ,  l2n  T  T  T  l  r, 

- f-  —  M  — ~  -Me  «J*  H  4-  peEx 

'be  3  ,:bC 


,  iy  < Pv 

Po—  =  P  — T  +  J*  Hx  +  Pe  E 
<'t 


p0—  -  ^  TT  +  ny  -  H*  +  P*E* 

' 't  ,  be 


It  is  assumed  that  the  nonlinear  perturbation  terms  are  negligible  in  comparison  with 
the  linear  terms.  Therefore,  the  energy  equation  becomes 

/)T'  /'u  -,2T' 

PoCv-^r  =  -  Po^o -r-  +  K—  (17) 

' 't  'be  i»x‘ 

Decoupled  Modes  of  Motion 

In  case  of  a  neutral  plasma,  i.e.,  pe  0,  the  number  of  ions  and  electrons  per 
volume  of  plasma  are  nearly  equal.  For  this  case,  if  one  considers  the  fact  that 

-  'hx  ' 'hx 

- =  0,  -  =  0,  then  it  is  possible  to  distinguish  between  two  modes  of  wave 

-  *t  i  be 

propagation,  transverse  mode  (z-direction)  and  longitudinal  mode.  In  the  transverse 
mode  the  states  are  found  to  be  ht  and  w,  and  the  state  equations  can  be  formed  from 
reduction  of  equations  (2), (3), (7), (8), (9), (l  1),  and  (16).  The  rest  of  the  equations  can  be 
reduced  to  form  state  equations  for  the  longitudinal  mode. 


:»S.» 


Transverse  Mode 

The  magneto-gas-dynamic  assumption  results  in  insignificant  magnetic 
induction  effect  in  Maxwell's  equations  from  the  terms  carrying  variations  of 
electric  field  with  time.  This  is  due  to  the  fact  that  nondimensional  parameters 
t0U  E0 

Rt  =  — —  and  Re  =  — —  are  of  the  order  of  one  or  smaller,  and 
L  MeLHo 

U2 

Rc  =  — ~  =  CJ2/iee  «  1  j  1 2 1 .  The  resulting  equations  are 


—  4-  H  — 

2  +  ■  i 

c  <  he 


V 1  -ffi 

V  X  uz 


where  i/H  =  -  and  V  =  /  — —  H.  The  parameter  V_  ‘is  defined  as  x- 

VPo 

component  of  speed  of  Alfven  wave. 


(ii.)  Longitudinal  Mode 

The  state  equation  for  this  mode  can  be  reduced  to 


H  .v 


i'(U  (’/y 

Hv  —  +  Hv  — 


'h 

x  lly 


=  -  RT, 


'£  f 

■he  -  *x 


4  «'2u  V2v  'h 

-f  —  v  — -  -  — 1 - - 


K  'i*T 


R  >  hi 

c„  1  *x 


»  n  ii  T  M e 

where  p  =  -c— ,  T  =  — ,  Vy  =  \  /  —  Hy.  The  parameter  V  is  defined 


y-component  of  the  speed  of  Alfven  wave. 


Lyapunov  Functional  and  Stability  Analysis 

In  this  section  the  Lyapunov  Functional  approach  is  applied  to  each  mode  of  the 
plasma  dynamics.  The  stability  results  are  derived  and  discussed. 


(i.)  For  the  transverse  mode  the  general  operator  form  of  the  equation  (18)  and  (19) 
evolution  equation  would  be 


where 


Z  -  Z(x.t) 


z= 


w 


'Z 


—  =  A(Z) 

< 't  ~ 


A  = 


/  r 


v2  /» 


Hx  ,8c 


(25) 


H. 


it 


bC 


.Assuming  that  the  boundary  conditions  are 

Z(O.t)  =0  ,  Z(f  ,t)  =  0 

The  operator  A  is  defined  in  L2(0/  )  and  its  domain  belongs  to  a  Hilbert  space.  A 
Lyapunov  functional  V(  Z(x,t),  t)  should  be  constructed  such  that 
dV 

V(Z,  t)  >  0  and  — —  <  0.  For  this  problem  the  functional  V  is  chosen  to  be  the 
dt 

inner  product  of  Z. 


V  =  <Z.Z>,  =  <Z,  PZ>0  ,  P  = 


or 


l  0 


0  or. 


,  or,,  or2  >  0. 


V  =  a,llhtll2  +  a2llwll2.  Hence;  V  =  —  <  Z,  Z  >,  =  2<  Z,  Z  >,  =  2<  AZ,  Z 

dt 


Since  P  is  symmetrix  <x.  Py>  =  <Px,  y> 

V  =  2<AZ,PZ>0  =  2<PAZ,Z>0 

where, 

r 

]2 


<  PAZ,  Z  >0  =  /  ZT  • 


,'ix! 


a,Hx  — 
1  'be 


V2  /, 
a - 

Hx  -be 


be2 


Z  dx 


(26) 


Considering  all  of  the  terms  in  the  above  integral,  the  following  can  be  resulted. 
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V=2  / 


■A 


aiL/n 


'He 


£  „  Hv 

+ 


He 


+  W 


«2v;  <'h 


£ 


Hv 


He 


+  Oc2i; 


w 


i  be2 


dx  (27) 


For  this  problem  since  the  eigenfunctions  of  states  are  similar,  i.e.,  the  Hilbert 
spaces  of  states  have  the  same  sequence  of  coordinates,  if  aj  andor2  are  selected 
such  that  =  V2/H 2  then  (27)  can  be  reduced  to 


V  =  2  / 


h 


'  ,2hT 


z  He2 


+  ot2v  w  ■ 


A 


dx  +  2alHxf 


i'Hv 


A 


h,~ - h  w  — 

-  be  rb c 


dx  (28) 


Using  integral  by  parts  and  the  following  relationships  [  1 3 j , 


/ 

J  1  He 


o  ( 

and 

c 


/ 


dx  >  tt2  J  f2  dx 
--'h 


,  '.Hv  ,  “£ 

h, -  +  w - 

z  'be  'be 


dx  =  hzw  j/  =  0 


the  following  can  be  resulted. 


V  =  -  2  / 


aluU 


A 

i be 


4-  <i7v 


.'Hv 
•  He 


dx 


<  —  27T2f(<yli/H  h2  +  0'2iAv2)dx 


if  «2  =  1,  etj  =  V2/H2  then 

|2 


Vx  2 

V  =  — y  Hhjll2 
H2 


+  llwl 


V  <  -2^ 


HJ 


Vft  Mh z 1 1  +  t/llwll2 


(29) 


In  this  case  (29)  shows  that  V  is  negative  definite  and  that  proves  the  stability  of 
the  transverse  mode  using  Lyapunov  approach. 


(ii.)  For  longitudinal  mode  of  wave  propagation  the  generic  form  of  evolution  equation 
(25)  is  considered 


r-'.y. 


i 


$ 


Z=  u  cL*(0/)  and  Z(0,t)  =  Z(f  ,t)  =  0 


In  this  case  A  in  equation  (25)  is  a  linear  operator  with  the  domain  in  a  separable 
Hilbert  space  of  state  function  which  maps  Jt  onto  itself.  From  equations  (20)  to 
(24)  A  is  formed  as. 


Vx2  <> 


Hx  <"‘x 


VV  / > 


Hy  he 


4 

-V— 7 
3  ,  he 


-RT„— 


r  ;> 


Cv  'he 


K  r)2 
Po^v  'he2 


An  approach  similar  to  the  transverse  mode  is  taken  to  construct  the  Lyaponov 
functional  for  the  longitudinal  mode. 


V  =  <Z,  PZ> 


where 


«2  O 

^3 


>  o 


Introduction  of  this  form  for  V  results  in  integral  terms  in  V  in  the  following  form 

/  Z,  —i  dx  . 


Such  terms  would  make  V  indefinite  unless  solution  forms  are  assumed  for  the  state 


variables.  In  order  to  make  the  stability  results  independent  of  the  solution  forms,  the 
following  equalities  are  assumed. 


v: 


aiHX  (-*2  TT 


H. 


-«lHy  =  -«3 


V 


,'2 


H„ 


-RT0ar3  =  — or4 


-RT0a3  - - or, 


R 


These  result  in  values  for  a2  through  c*4  and  force  all  of  the  integral  terms  in  form 
of  (30)  to  go  to  zero.  Assumine  oq  is  unity  the  following  are  resulted. 


A 


rtj  =1 


H, 


«2  - 


,-2 


a3  = 


V 


h; 


«4  =  77^  RT 

V 


H‘ 


^  =  y7CVT 


Based  on  these  values  for  <*,’s  and  similar  algebraic  techniques  used  in  the  transverse 
mode,  the  following  can  be  derived. 


V  =  V  or,  I IZ,  1 12  >0 

i- 1 

and 

V  <.  —'In2  (*,//Hllhyll2  +  o2z^l  Iv  1 12  4-  r»3—  i/llu|  |2  4-  a4 — — — I  IT"  1 12 

3  P  o^v 

These  results  indicate  that  for  Z,  &  nr,  *  0  V  is  positive  and  V'  is  negative. 
Therefore,  the  longitudinal  mode  is  stable.  Hence,  stability  of  the  MPD  engine  near  its 
equilibrium  is  established  without  solving  any  of  its  dynamic  equations. 

Stability  results  only  for  the  transverse  mode  of  the  MPD  engine  can  be  established 
by  application  of  point  spectrum  approach.  This  would  provide  a  means  of  checking 
the  results  presented  in  this  paper,  namely,  by  means  of  the  Lyapunov  functional 
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approach. 

By  means  of  separation  of  variables  the  semi-group  property  T(t),  generated  by  A 
can  be  constructed. 


or 


T, 


H, 


Xh 


z  = 


-=r  -  ^HTT-  H, 

1  h  Ah 


ThXh 

TWXW 


Xv 

x„ 


T 

Tu 


u 


Xu 


ThXh 

TWXW 


=  0 


K 

,-v 


i 


If: 


The  solution  for  Z  = 
i.e., 


ThXh 
T  X 

1  wvw 


will  exist  if  the  above  operation  on  Z  is  not  one-to-one; 


det 


\  Xh 

h  Ah 

XL 

Hv  Xh 


Hv 


X 


w 


T  X 

i  w  ^vw 

-  v 


T 


w 


X 


w 


=  0 


(31) 


In  order  to  have  T  and  X  functions  independent  from  x  and  t,  respectively,  it  is 

T  X"  T  x'  X 

required  that  —  =*=  f(x,t),  ~  =*  f(x,t).  From  — ,  77  =  const,  and—  =  const,  it  is 
T  X  T  X  X 

conclusive  to  represent  X  functions  in  terms  of  a  real  periodic  function 


X  =  Are,Xx/'  +  A, 


-  iXx/f 


where  A  =  complex  conjugate  of  A 


X_ 

X 


f2 


x_ 

X 


iX 

(~ 


and  T  =  Afe31  — *  —  =  S.  From  the  boundary  condition  it  appea 


rs 


X  =  Ax  sin  — 


£1 


Therefore,  equation  (31)  can  be  reduced  to 


det 


xn2 

XL  i 

H,  “  f 


H  i  — 
Hx 1  y 


S  +  v 


( 


0 


The  point  spectrum  of  operator  A,  i.e.,  crp(A)  can  be  found  as 
Op(A)  =  jS  f  of  A)  [  (SI— A)  is  not  one  to  one  ) 


S2  +  S(i/H  +  v)  ~J~  +  Vx2  yy  +  i/V„  —■  =  0 


The  roots  of  this  equation  are 


I  1,2 


1 

2 


-  {"H  +  ») 


( 


V 


(^H  - 


e 4 


-  4V‘ 


X2 

An 

2 


The  above  expression  for  Sn  indicates  that  sup  [Re  a  (A))  <  0  which  is  the  necessary 
and  sufficient  condition  for  the  equilibrium  solution  of  equations  (18)  and  (19)  to  be 
exponentially  stable,  i.e.  an  equivalence  to  uniform  assymptotic  stability  for  the  above 
linear  system. 


The  point  spectrum  approach  applied  to  the  longitudinal  mode  results  in  the 
following  spectrum  of  the  set  of  state  equations. 


K(—  +  {  —  S)X4  + 
Po  3  p0 


S2-L  +  1  t/S2 

Pc  3  T0(7-l) 


+  scp 


s3 

T0(7-l) 


(^HX2  +  S)(NuX2  +  S)  +  V2X2 


+  X2V2(S  +  iA2) 


TQ(7-l)i 


S  K  X2 
Po 


0 


This  characteristic  equation  doesn't  have  a  closed  form  solution.  Therefore,  in 
general  the  spectrum  approach  doesn’t  yield  stability  solution  for  a  system,  unless  in 
very  simplified  and  special  cases.  Whereas  the  outlined  Lyapunov  approach  would 
result  in  stability  solutions  for  a  system. 


Conclusion 


A  new  approach  to  stability  analysis  of  an  MPD  engine  was  presented.  This 
technique  is  based  on  the  Lyapunov  stability  analysis  which  is  extended  to  cover 
distributed  parameter  systems.  The  results  of  the  stability  analysis  was  supported 
with  those  derived  from  specteral  analysis  point  of  view.  The  procedure  for 
construction  of  the  Lyapunov  functional  and  derivation  of  its  derivative  was 
presented.  It  was  shown  that  while  spectral  technique  can  be  applied  to  stability 
analysis  of  a  special  class  of  systems,  the  presented  method  can  be  applied  to  any  form 
of  distributed  parameter  systems. 

Nomenclature 


Constant  volume  specific  heat 
Const,  press,  specific  heat 
Internal  energy/mass 
Electric  field  vector 

Magnetic  field  vector,  perturbation  vector  of  magnetic  field. 
Constant  magnetic  field  or  characteristics  magnetic  field 
Electric  current  density 
Thermal  conductivity  of  plasma 

Characteristic  length,  length  in  the  x  direction  (for  the  flow) 
Pressure,  perturbation  in  pressure 
Heat  diffusion/area 
Gas  constant 

Temperature,  perturbation  in  temperature 

time,  characteristic  time 

Velocity  vector,  characteristic  velocity 

velocity  vector  components  in  x,y  and  z  directions 

Spatial  coordinates 

Spatial  coordinates  for  j  =  1,  2,  3 


P,  p: 
Q: 


u,  v,  w: 


x,  y,  z: 


Greeks 


Dielectric  constant 
Viscosity,  kinematic  viscosity 
Magnetic  permeability 
Plasma  density,  perturbation  in  density 


V.V.V.V'^' 
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pe:  Excess  electric  charge/volume 


cr. 


uH: 


Electrical  conductivity 

1 

ane 


subscripts 

x,  y,  z:  Variable  subscribed  in  the  direction  of  x,  y  and  z. 

Notation 

Th,  Tw:  Time  dependent  part  of  variable  in  subscript 
X^,  Xw:  x  dependent  part  of  variable  in  subscript. 
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